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ABSTRACT 

Unmanned  Air  Vehicles  have  become  increasingly  important  on  the  modern 
battlefield.  The  restrictive  requirement  for  runways  and  special  equipment  to  take  off 
and  land  was  partially  solved  by  the  vertical  take  off  and  landing  Airborne  Remotely 
Operated  Device,  AROD.  Work  done  at  the  Naval  Postgraduate  School  has  modified 
the  AROD  to  not  only  land  and  launch  vertically,  but  to  fly  horizontally  for  the 
majority  of  the  mission.  To  realize  these  capabilities,  as  well  as  that  of  autonomous 
flight,  an  accurate  computer  model  was  required  of  both  the  AROD  and  the  avionics 
test  bed  aircraft.  Bluebird,  in  order  to  design  the  control  and  navigation  systems.  High 
fidelity,  non-linear  equations  of  motion  were  derived  in  matrix  form  that  represented 
any  six  degree  of  freedom  aircraft  model,  and  were  then  tailored  for  use  on  specific 
aircraft.  Computer  modeling  of  the  resulting  equations  of  motion,  as  well  as  the 
sensors  used  on  the  aircraft,  was  done  using  SiMULINK  and  Matlab  software.  The 
resulting  computer  model  provides  a  non-linear  system  of  equations,  which  are  easily 
linearized  at  any  desired  flight  condition,  as  required  by  the  proposed  control  and 
navigation  system  design. 


m 


TABLE  OF  CONTENTS 

I.  INTRODUCTION 1 

A.  IMPORTANCE  OF  UNMANNED  AIR  VEHICLES 1 

B.  UAV  RESEARCH  USING  THE  AROD  VEHICLE 2 

C.  REQUIREMENT  FOR  MODELING 3 

II.  BACKGROUND 4 

A.  DESCRIPTION  OF  AROD 4 

B.  AROD  MODELING 5 

C.  DESCRIPTION  OF  BLUEBIRD 8 

III.  AIRCRAFT  EQUATIONS  OF  MOTION 9 

A.  NOTATION 9 

B.  COORDINATE  SYSTEMS 10 

C.  SPATIAL  ORIENTATION 11 

1.  Euler  Angles 11 

2.  Quaternions 15 

D.  DERIVATION  OF  EQUATIONS  OF  MOTION 18 

1.  Linear  Equations 19 

2.  Angular  Equations 20 

3.  State  Equations 21 

E.  EXTERNAL  FORCES  AND  MOMENTS 22 

1.  Aerodynamic  Forces  and  Moments 22 

2.  Other  Forces  and  Moments 25 

IV.  COMPUTER  MODELING  OF  THE  AIRCRAFT  EQUATIONS  OF 
MOTION 28 

iv 


A.     IMPLEMENTATION  OF  EQUATIONS 


WNTEREVci^-pSC 


scHoo; 


1.  Procedure 28 

2.  AROD  Equations 29 

a.  Control  Forces  and  Moments 30 

b.  Additional  Forces  and  Moments     31 

3.  Model  Validation 32 

a.  Kinematic  Equations     32 

b.  Gravitational  Forces 34 

c.  Additional  Forces  and  Moments     36 

4.  Bluebird  Equations  of  Motion 43 

a.  Kinematic  Equations     43 

b.  Gravitational  Forces 45 

c.  Aerodynamic  Forces  and  Moments 47 

B.     VALIDATION  OF  AN  INDEPENDENT  CASE 50 

V.        SENSOR  AND  ACTUATOR  MODELING 56 

A.  ACCELEROMETER  MODELING 56 

1.  Error  Model     59 

2.  Results  and  Validation 61 

B.  RATE  GYRO  MODELING 64 

1.  Error  Modeling 65 

2.  Results  and  Validation 66 

C.  PITCH,  ROLL,  AND  HEADING  SENSOR  MODELING    70 

1.       Error  Modeling 70 

D.  MODELING  OF  AN  ARBITRARY  SENSOR  PLACEMENT   ...  77 

1.  Linear  Accelerations 77 

2.  Angular  Accelerations 79 

V 


\ 


VI.       CONCLUSIONS  AND  RECOMMENDATIONS 80 

A.  CONCLUSIONS 80 

B.  RECOMMENDATIONS 80 

APPENDIX  A:         MATHEMATICAL  PROPERTIES 82 

A.  CROSS  PRODUCT  PROPERTIES    82 

B.  DERIVATIVES  OF  VECTORS     83 

C.  EQUATIONS  USED  FOR  LINEARIZATION 84 

APPENDIX  B:           NUMERICAL  RESULTS 86 

A.  AROD  RESULTS 86 

B.  BLUEBIRD  NUMERICAL  RESULTS 90 

C.  CESSNA  172  RESULTS 94 

APPENDIX  C:         PROGRAM  LISTINGS 96 

A.  AROD  Matlab  ROUTINES 96 

1.  Main  Routine 96 

2.  Supporting  Subroutines     100 

3.  Data  and  Initialization  Subroutines 102 

B.  BLUEBIRD  Matlab  ROUTINES 106 

1.  Main  Routine 106 

2.  Supporting  Subroutines     113 

3.  Data  and  Initialization  Subroutines 115 

REFERENCES 120 

INITIAL  DISTRIBUTION  LIST 122 


VI 


LIST  OF  TABLES 

2.1  PHYSICAL  CHARACTERISTICS  OF  AROD 6 

2.2  VANE  DEFLECTION  COMBINATIONS  FOR  POSITIVE  ANGLES  6 

2.3  PHYSICAL  CHARACTERISTICS  OF  Bluebird 8 

4.1  NON-DIMENSIONAL  DERIVATIVE  DATA  FOR  AROD 42 

4.2  NON-DIMENSIONAL  STABILITY  DERIVATIVES 49 

4.3  NON-DIMENSIONAL  CONTROL  DERIVATIVES 50 

4.4  COMPARISON  OF  EIGENVALUES  FOR  BLUEBIRD 50 

4.5  EIGENVALUE  COMPARISON  FOR  CESSNA  172  TEST  CASE    .  .  51 

5.1  ACCELEROMETER  CHARACTERISTICS 60 

5.2  GYRO  CHARACTERISTICS 66 

5.3  INCLINOMETER  AND  HEADING  SENSOR  CHARACTERISTICS  71 


Vll 


LIST  OF  FIGURES 

2.1  Airborne  Remotely  Operated  Device 5 

2.2  AROD  Direction  of  Positive  Vane  Deflections     7 

3.1  Relative  Position  of  Coordinate  Systems 9 

3.2  Z-Y-X  Euler  Angle  Rotation  Sequence 12 

4.1  Gyroscopic  Motion  of  AROD     33 

4.2  Modeling  of  Gravitational  Effects  for  AROD 34 

4.3  Velocity  of  AROD,  Gravitational  Effects  Included 37 

4.4  Thrust  vs.  RPM  for  AROD 38 

4.5  Thrust  and  Moment  Data  for  AROD 38 

4.6  Moment  Due  to  Vane  Deflection 40 

4.7  Non-Linear  Equations  of  Motion  Model 42 

4.8  Rolling  Motion  For  Complete  Non-Linear  Equations     43 

4.9  Block  Diagram  of  Kinematic  Equations  of  Motion 44 

4.10  Gravitational  Forces  Model     46 

4.11  Gravitational  Effects  on  Velocity     48 

4.12  Full  Non-Linear  Equations  of  Motion  Model 49 

4.13  Comparison  of  Longitudinal  Responses  to  Step  Elevator  Input    ....  52 

4.14  Comparison  of  Lateral  Responses  to  Step  Rudder  Input 53 

4.15  Comparison  of  Lateral  Responses  to  Step  Aileron  Input 53 

4.16  Difference  in  Analytic  and  Numerical  Results,  Step  Elevator  Input  .  .  54 

4.17  Difference  in  Analytic  and  Numerical  Results,  Step  Rudder  Input     .  .  54 

4.18  Difference  in  Analytic  and  Numerical  Results,  Step  Aileron  Input     .  .  55 
5.1  Typical  Accelerometer  Model     57 

viii 


5.2  Accelerometer  Modeling 58 

5.3  Synthesis  Model  for  Accelerometers 59 

5.4  Error  Model  for  Accelerometers 61 

5.5  Measured  and  True  Acceleration  From  a  Step  Elevator  Input 62 

5.6  Measured  and  True  Acceleration  From  a  Step  Aileron  Input 63 

5.7  Measured  and  True  Acceleration  From  a  Step  Rudder  Input 63 

5.8  Functional  Diagram  of  a  Rate  Gyro 65 

5.9  Rate  Gyro  Model 66 

5.10  Rate  Gyro  Synthesis  Model     67 

5.11  Measured  and  True  Angular  Rates  From  a  Step  Aileron  Input     ....  68 

5.12  Measured  and  True  Angular  Rates  From  a  Step  Elevator  Input  ....  68 

5.13  Measured  and  True  Angular  Rates  From  a  Step  Rudder  Input    ....  69 

5.14  Simple  Pendulum  Inclinometer 71 

5.15  Inclinometer  Error  Modeling 72 

5.16  Response  to  a  Step  Elevator  Input 75 

5.17  Response  to  a  Step  Rudder  Input 75 

5.18  Response  to  Step  Aileron  Input 76 


IX 


ACKNOWLEDGMENT 

I  would  like  to  express  my  appreciation  for  the  assistance  given  in  the  develop- 
ment of  this  thesis  by  Dr  I.  I.  Kaminer  and  by  Dr.  R.  M.  Howard.  Their  boundless 
patience  and  professional  counsel  was  invaluable  to  the  sucessful  completion  of  this 
project. 


I.  INTRODUCTION 

A.     IMPORTANCE  OF  UNMANNED  AIR  VEHICLES 

Unmanned  Aerial  Vehicles  have  become  increasingly  more  important,  both  on 
the  battlefield  and  in  civilian  service,  since  the  Ryan  Q-2C  Firebee  target  drone 
introduced  the  "modern  age"  of  UAVs  in  1960  [Ref.  Siu  91].  From  that  time  on, 
military  planners  have  assured  that  UAVs  have  the  capability  to  collect  intelligence, 
target  enemy  positions,  gather  bomb  damage  assessment,  as  well  as  perform  many 
other  tasks.  The  real  benefit  in  using  unmanned  aircraft  lies  in  the  fact  that  many 
missions  can  be  performed  deep  in  enemy  territory,  all  without  endangering  the  lives 
of  pilots,  or  risking  the  loss  of  a  much  more  expensive  aircraft.  With  the  recent  use  of 
UAVs  in  Operation  Desert  Storm,  improvements  in  the  current  technology  are  both 
indicated  and  desirable. 

The  most  capable  UAV  in  service  of  the  United  States  Navy  and  Marine  Corps 
today  is  the  Pioneer  Short  Range  UAV.  The  system  is  hampered,  though,  by  the 
large  amount  of  runway  and  special  equipment  needed  to  launch  and  land  the  aircraft. 
These  requirements  limit  the  usefulness  of  the  Pioneer  by  keeping  the  aircraft  take- 
off and  landing  area  well  away  from  the  areas  where  the  ground  forces  are  operating. 
This  distance  then  leads  to  longer  transit  times  to  and  from  the  assigned  operating 
area,  and  thus  a  shorter  time  on  station.  However,  what  is  really  required  in  many 
instances  by  the  ground  forces  is  an  aircraft  that  can  respond  quickly  to  a  changing 
tactical  situation.  The  Airborne  Remotely  Operated  Device,  AROD,  was  an  attempt 
by  Sandia  National  Laboratories  [Ref.  Wh  87],  in  response  to  a  requirement  by  the 
Naval  Ocean  Systems  Center,  NOSC,  to  respond  to  these  needs. 


The  United  States  Marine  Corps  had  set  a  requirement  for  a  short  range,  direct 
support  UAV  as  described  in  [Ref.  MCG  87]: 

"  . . .  to  allow  the  front  line  commander  to  see  "over  the  next  hilF",  to  a 
distance  of  two  kilometers  ..." 

The  AROD  was  designed  to  be  a  ducted  fan,  hovering  device  carrying  a  fiber  optic 
data  link  and  on  board  cameras.  AROD  testing  was  canceled  [Ref.  Sa  89]  as  the 
Department  of  Defense  requirements  grew,  requiring  a  minimum  range  of  30  km  for 
all  Short  Range  UAVs.  The  AROD  was  incapable  of  this  kind  of  range,  however;  the 
design  is  still  potentially  useful. 

The  Unmanned  Air  Vehicle  Flight  Research  Lab.  UAV  FRL,  at  the  Naval  Post- 
graduate School  has  proposed  a  solution  using  the  AROD  that  would  satisfy  the 
DoD  short  range  UAV  requirements,  while  maintaining  the  important  capability  for 
vertical  take-offs  and  landings. 

B.     UAV  RESEARCH  USING  THE  AROD  VEHICLE 

Unmanned  Aerial  Vehicle  research  underway  at  NFS  has  taken  the  AROD  air- 
frame and  fitted  it  with  wings  from  the  .^quila  UAV  [Ref.  Kre  92,  Sto  93]  in  order 
to  give  the  AROD  forward  flight  capability.  The  proposed  configuration  will  give  the 
Archytas  (an  AROD  with  wings)  the  ability  to  take-off  and  land  vertically  and  then 
transition  to  horizontal  flight  for  the  mission.  This  design  will  explore  new  technol- 
ogy, driven  by  the  goals  established  by  the  UAV  Joint  Project  Office  [Ref.  DOD  92] 
of: 

•  Take  off  weight  under  200  lb 

•  Carry  a  50  lb  payload 

•  Fly  at  a  maximum  speed  of  150  kts 


•  Take-ofF  and  land  in  an  area  30m  by  60m 

The  vertical  flight  is  accomplished  with  a  powerful  ducted  fan, which  causes  a  great 
deal  of  gyroscopic  coupling  and  torque  when  producing  enough  thrust  to  lift  the 
aircraft.  Therefore,  in  order  to  achieve  stable  take-offs  and  landings,  a  three-axis 
autopilot  is  a  necessary  feature  of  the  aircraft.  Additional  capabilities  desired  in 
the  final  version  of  the  Archytas  are  guidance  and  navigation  systems  which  will 
allow  autonomous  operation,  as  well  as  a  Global  Positioning  System  aided  autoland 
capability. 

C.     REQUIREMENT  FOR  MODELING 

Simulation  and  modeling  of  the  aircraft  are  essential  to  the  successful  design 
of  a  control  system  capable  of  autonomous  flight.  The  model  must  be  a  very  high 
fidelity,  non-linear  model,  that  can  be  easily  linearized  at  any  given  flight  condition. 
The  model  should  be  able  to  interpolate  between  data  points  resulting  from  wind 
tunnel  testing  in  order  to  simulate  the  highly  non-linear  transition  from  vertical  to 
horizontal  flight.  Moreover,  the  model  must  also  be  capable  of  including  the  outputs 
of  the  sensors  as  inputs  to  the  control  and  navigation  system  for  sensors  located  at 
any  arbitrary  location  on  the  aircraft. 

This  thesis  develops  a  six  degree  of  freedom  model  for  the  A  ROD  in  the  vertical 
flight  regime,  as  well  as  for  an  aircraft  in  a  fixed  wing  configuration.  This  test  aircraft, 
named  Bluebird,  is  used  to  test  the  guidance,  navigation,  and  control,  GNC,  systems 
in  horizontal  flight,  since  there  currently  is  no  aerodynamic  data  available  for  the 
Archytas  configuration.  Use  of  the  Bluebird  will  provide  the  capability  to  design  and 
test  a  GNC  system  on  a  stable  aircraft  before  the  first  flight  on  the  Archytas. 


II.  BACKGROUND 

A.     DESCRIPTION  OF  AROD 

The  AROD  was  designed  by  the  Sandia  Research  Laboratory  in  Albuquerque, 
New  Mexico  in  a  project  managed  by  NOSC.  The  vehicle  possessed  no  flying  surfaces 
and  relied  solely  on  powered  lift  for  flight.  Control  of  the  aircraft  was  obtained  through 
the  use  of  four  fixed  anti-torque  vanes  and  four  moveable  control  vanes  positioned  in 
the  propeller  wash  of  the  duct  [Ref.  We  88].  The  main  features  of  the  AROD  were 
vertical  take-off  and  landing,  VTOL,  flight,  lightweight  construction,  compact  size, 
and  minimal  support  equipment  required.  However,  the  AROD  required  most  of  the 
engine  output  to  maintain  the  powered  lift,  so  very  little  excess  thrust  was  left  for 
translational  flight. 

An  important  aspect  of  the  AROD  design  was  the  improvement  in  static  perfor- 
mance provided  by  the  efficiency  of  the  ducted  fan  design.  The  addition  of  the  shroud 
around  the  three-bladed  propeller  resulted  in  increased  mass  flow  through  the  fan, 
and  thus  more  static  thrust  when  compared  to  a  conventional  propeller  configuration 
[Ref.  Kre  92].  The  AROD  is  shown  in  Figure  2.1  and  characteristics  of  the  AROD 
are  tabulated  in  Table  2.1.  The  moveable  control  vanes  are  all  used  in  combination 
to  exert  the  desired  control  forces  on  AROD.  Roll  control  is  achieved  by  deflecting  the 
four  vanes  in  the  same  direction,  while  pitch  and  yaw  control  is  obtained  by  deflecting 
a  pair  of  vanes  in  the  required  direction.  The  numbering  of  the  vanes  is  shown  in 
Figure  2.2  and  the  combinations  of  vane  deflections  required  for  positive  roll,  pitch, 
and  yaw  motion  are  given  in  Table  2.2. 
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Figure  2.1:  Airborne  Remotely  Operated  Device,  [Ref.  Siu  91] 

B.     AROD  MODELING 

The  AROD  vehicle  has  been  the  subject  of  several  theses  at  NPS.  The  designs 
presented  in  the  theses  rely  on  the  AROD  model  given  by  Sandia  Labs  in  the  original 
design  of  their  controller  [Ref.  Wh  87,  Wh  91].  This  model  was  based  on  the  more 
classical  technique  of  linearizing  an  aircraft  model,  based  on  dimensional  derivatives 
in  a  state  space  form.  The  resulting  model  was  acceptable  for  the  AROD  in  a  hovering 
and  near  vertical  translational  flight  mode,  but  was  not  easily  adaptable  to  anything 
other  than  the  narrow  range  of  conditions  planned  for  AROD.  The  Sandia  Labs 
papers  also  pointed  out  several  types  of  coupling  in  the  AROD.  The  most  prominent 
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TABLE  2.1:  PHYSICAL  CHARACTERISTICS  OF  AROD 


Inlet  Diameter,  A 

29.25  in 

Propeller  Radius,  R 

12  in 

Exit  Radius 

23.375  in 

Inlet  Area  Ratio 

1.219 

Exit  Area  Ratio 

1.115 

Exterior  Contour 

Tapered  Rear 

Propeller  Location,  %  chord 

25% 

Number  of  Blades 

3 

Engine  Speed,  Max. 

8000  rpm 

Engine  Speed,  Nom. 

6500  rpm 

Tip  Speed,  Max. 

838  fpm 

Tip  Speed,  Nom. 

680  fpm 

Power  Loading,  ^^g^^"' 

7.25  HP/p 

Mass  Moment  of  Inertia,  /j- 

1.2312  slug-  p 

Mass  Moment  of  Inertia.  ly 

3.9548  slug  -  p 

Mass  Moment  of  Inertia,  /^ 

3.9825  slug  -  p 

Prop  Mass  Moment  of  Inertia,  In 

0.00898  slug  -  p 

Prop  Mass  Moment  of  Inertia,  Iry 

0.0045  slug  -  p 

Prop  Mass  Moment  of  Inertia,  7^ 

0.0045  slug  -  p 

TABLE  2.2:    VANE  DEFLECTION  COMBINATIONS  FOR  POSITIVE 
ANGLES 


Vane  Combination 

Roll,  $ 

Vi  +  V2  +  V3  +  V4 

Pitch,  0 

V  2  -  V4 

\'a\v,  vl/ 

Vi  -  V3 

of  the  coupling  effects  is  the  gyroscopic  coupling  between  the  pitch  and  yaw  axes 
resulting  from  the  large  amount  of  angular  momentum  contributed  to  the  aircraft 
by  the  propeller.  Another  dynamic  coupling  exists  between  the  altitude-rate  and  the 
vehicle  attitude,  since  a  loss  of  lift  due  to  thrust  will  occur  when  the  vehicle  is  tilted  to 
generate  horizontal  motion,  ^'et  a  third  dynamic  coupling  exists  between  the  altitude 


Figure  2.2:  AROD  Direction  of  Positive  Vane  Deflections 


and  roll  control  loops,  since  the  reactive  torques  applied  to  the  roll  axis  vary  as  the 
engine  speed  is  varied.  Sandia  Labs  also  provided  data  for  modeling  both  the  engine 
and  the  servos  as  second  order  transfer  functions  which  were  used  in  this  thesis. 

Additional  information  was  obtained  by  Weir  [Ref.  We  88]  in  wind  tunnel  test- 
ing. This  information  included  non-dimensional  derivatives  for  vane  effectiveness  and 
non-dimensional  stability  derivatives.  The  report  also  stated  that  the  control-vane 
effectiveness  is  constant  out  to  at  least  25  deg  of  deflection.  Wind  tunnel  data  were 
also  presented  to  show  that  control  vane  effectiveness  in  approximately  the  same  for 
translational  flight  as  for  hovering  flight.  This  equivalence  is  due  to  the  fact  that  the 
vanes  are  located  in  the  high  speed  flow  aft  of  the  propeller  and  are  not  significantly 
affected  by  the  freestream. 


TABLE  2.3:  PHYSICAL  CHARACTERISTICS  OF  BLUEBIRD 


Weight 

55  lbs 

Average  Wing  Chord,  c 

1.802  / 

Wing  Span,  b 

12.42  / 

Planform  Surface  Area,  S 

22.380  p 

Engine  Power 

4.0  HP 

Mass  Moment  of  Inertia,  /j. 

10.0  slug-  p 

Mass  Moment  of  Inertia,  ly 

16.12  slug-  p 

Mass  Moment  of  Inertia,  I^ 

7.97  5/1/^  -  p 

C.     DESCRIPTION  OF  BLUEBIRD 

The  Bluebird  aircraft  was  acquired  as  a  test  bed  for  guidance  and  navigation 
systems.  Ultimately,  these  systems  will  be  installed  on  the  Archytas  aircraft.  The 
Bluebird  is  a  conventional  aircraft  that  will  be  used  to  test  systems  for  a  similar 
configuration  to  the  Archytas  in  forward  flight.  The  aircraft  model  was  developed 
in  the  same  manner  as  for  AROD,  as  will  be  described  in  Chapter  III.  Physical 
characteristics  for  the  Bluebird  are  given  in  Table  2.3. 


III.  AIRCRAFT  EQUATIONS  OF  MOTION 

A.     NOTATION 

In  this  section  the  notation  used  in  modeling  the  equations  of  motion  is  intro- 
duced. This  notation  is  common  in  the  field  of  robotics  (see  [Ref.  Sil  91]  and  [Ref. 
Cra  86]).  and  is  shown  below  in  Figure  3.1.  The  following  conventions  are  used: 


{A) 


<. 


^    "  BO 


Figure  3.1:  Relative  Position  of  Coordinate  Systems 


{y4}  represents  the  coordinate  system  with  basis  vectors,  x^,?/^,  and  za- 

^Pq  represents  the  position  of  point  Q,  expressed  in  {>4}. 

^Vq  represents  the  velocity  of  point  Q,  measured  in  {yl}  and  expressed  in  {A}. 

^{^Vq)  represents  the  velocity  of  point  Q,  measured  in  {A},  and  expressed  in 
{B). 


•  qR  is  a  rotation  matrix,  also  called  a  direction  cosine  matrix.    A  free  vector 
in  one  coordinate  system,  that  is  a  vector  that  can  be  moved  parallel  to  itself 


without  change  e.g.,  ^Vq  can  be  expressed  in  another  coordinate  system  by 
using  the  rotation  matrix: 

•  "^Ub  is  the  angular  velocity  of  the  {B]  coordinate  system  with  respect  to  {A}^ 
and  expressed  in  {A}. 

•  ^{^^b)  is  the  angular  velocity  of  {B].  with  respect  to  {A],  and  expressed  in 
{B}. 

•  Additional  information  can  be  added  to  the  subscripts  i.e.,  "^Pbo  is  the  position 
of  the  origin  of  [B],  expressed  in  {^4}. 

B.     COORDINATE  SYSTEMS 

In  order  to  derive  equations  of  motion  for  a  rigid  airplane,  the  following  coor- 
dinate systems  and  assumptions  will  be  used: 

•  [U]  represents  the  inertial  tangent  plane  coordinate  system  attached  to  Earth. 

•  {5}  represents  the  body  fixed  coordinate  system. 

•  All  sensors  are  located  at  the  e.g.  (This  assumption  will  be  lifted 
in  a  later  section) 

•  The  mass  of  the  aircraft  remains  constant. 

•  Given  a  vector  u,  its  derivative  with  respect  to  [B]  is  denoted  as  ^(•) 
and 

its  derivative  with  respect  to  [V]  is  denoted  as  (') 
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The  {B}  coordinate  system  is  defined  with  Xb  as  the  thrust  axis.  A  positive  roll 
rate,  p,  is  clockwise  when  looking  in  the  positive  X  direction.  The  positive  direction 
for  Yb,  the  pitch  axis,  is  out  the  right  wing  .  A  positive  pitch  rate,  q,  is  defined  as 
clockwise  looking  in  the  positive  Y  direction.  The  Zb  axis  is  the  yaw  axis,  and  a 
positive  yaw  rate,  r,  is  defined  as  clockwise,  looking  in  the  positive  Z  direction. 

To  simplify  the  notation  in  places  where  it  becomes  cumbersome,  The  following 
definitions  are  introduced: 

•  vq  represents  the  velocity  of  an  arbitrary  point,  Q,  measured  and  expressed  in 
{U]. 

•  vbo  represents  the  velocity  of  the  origin  of  {B},  measured  and  expressed  in 
{t/},  i.e.  ^'Vbo   =  VBO- 

•  ^vq  represents  the  velocity  of  point  Q,  measured  in  {U}  and  expressed  in  {5}, 
i.e.  ^{"Vq)   =  ^VQ. 

•  u)B  represents  the  angular  velocity  of  {B},  measured  and  expressed  in  {U],  i.e. 
^'Qb    =  i^B- 

•  ^ujb  represents  the  angular  velocity  of  {B}  measured  in  {U]^  and  expressed  in 
{B}/i.e.^{''nB)   =  ^uJB. 

C.     SPATIAL  ORIENTATION 
1.      Euler  Angles 

The  spatial  orientation  of  a  rigid  body  [Ref.  Ju  92]  can  be  defined  by  the 
three  Euler  angles,  ^,0,  and  ^  called  roll,  pitch  and  yaw  and  defined  in  Figure  3.2. 
The  Euler  angles, in  turn,  can  be  used  to  define  a  rotation  between  two  coordinate 
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Y   =Y 


Figure  3.2:  Z-Y-X  Euler  Angle  Rotation  Sequence 

systems.  This  rotation  is  obtained  using  Euler's  theorem: 

Any  number  of  rotations  about  different  axes  through  a  point  must,  in 
the  end,  remain  equivalent  to  a  single  rotation. 

For  the  case  of  conventional  aircraft,  a  3-2-1  rotation  sequence  is  used  [Ref.  Sch  92], 
where  the  aircraft  is  yawed,  pitched,  then  rolled.  In  the  cases  investigated  here,  0  is 
small,  and  in  steady  state  flight  is  equal  to  the  angle  of  attack,  q.  $  can  be  expected 
to  be  anywhere  from  ±  60deg  in  normal  flight  and  can  be  anywhere  from  ±  ISOdeg  in 
acrobatic  flight.  ^  represents  the  heading  of  the  aircraft  and  of  course  can  range  from 
0  to  360  deg.  The  transformation  from  inertial  coordinates{(/},  to  body  coordinates 
{5},  is  carried  out  as  follows,  and  is  shown  in  Figure  3.2. 

1.  The  inertial  coordinate  system  is  represented  by  the  vector  ^'\\  with  the  com- 
ponents X,  y,  and  z.  The  first  rotation  is  made  about  the  z  axis  through  an  angle 
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^.  Now  the  vector  is  expressed  as  ^V  with  the  components  I2,  y2^  and  Z2-  Since 
the  rotation  was  about  the  z  axis,  the  22  component  remains  unchanged.  The 
resulting  elemental  matrix  is: 


M{<lf)  = 


cos  ^      sin  ^P     0 

-  sin  ^    cos  ^    0 

0  0        1 


(3.1) 


2.  Now  the  rotation  is  made  about  the  new  y  axis,  1/2,  through  an  angle  0.  This 
results  in  a  third  coordinate  system  with  the  vector  expressed  as  ^V,  and  having 
components  x^^y^,  and  ^3.  This  rotation  does  not  change  the  ^3  component. 
The  resulting  elemental  matrix  is: 


M(0)- 


cos  0    0     —  sin  0 

0        1  0 

sin  0     0      cos  0 


(3.2) 


3.  Lastly,  the  rotation  about  the  2:3  axis  through  an  angle  0  is  made  to  give  the 
vector  expressed  in  body  coordinates,  ^V.  Now  the  resulting  matrix  is 


M{^)  = 


1  0  0 

0      cos  $      sin  $ 
0     -  sin  0    cos  $ 


(3.3) 


When  the  matrices  are  multiplied  together  in  the  correct  sequence,  M($)A/(0)A/(^), 

the  result  is  the  yR  direction  cosine  matrix,  expressed  in  terms  of  Euler  angles  as 

shown 

cos^  COS0  sin^  cos0  —  sin0 


cos^  sin0sin$  —  sin^  cos<E>    sin0  sin$  sin^  +  cos^  cos$     cos0sin$ 
cos^  sin0  cos$ -I- sin^  sin$    sin0cos<^sin^  —  cos^  sinO    cos0cos$ 


(3.4) 


The  next  step  is  to  develop  the  kinematic  differential  equations  that  de- 
scribe the  change  in  Euler  angles  with  time.  Following  the  method  used  in  [Ref. 
Sch  92],  the  matrix  of  differential  equations,  J7,  can  be  written  as  a  sum  of  individual 
Euler  angle  rates: 


n  =  M($) 


0 

0 

+  M(O)M(0) 

i^ 

0 

0 

+  M(4>)M(0)M(^) 


0 
0 


(3.5) 


13 


When  the  matrix  algebra  in  Equation  3.5  is  done,  the  resulting  kinematic  differential 
equations  for  p,^,and  r  are  given  as: 


P 
Q 

r 


1  0  -sinG 

0      cos  ^       cos  0  sin  $ 
0    —  sin  $    cos  0  cos  $ 


0 


(3.6) 


The  matrix  on  the  right  hand  side  of  Equation  3.6  is  invertible  for  all  0  ^  7r/2,  and 
can  be  used  to  solve  for  the  Euler  angle  rates,  $,  0  and  ^: 


0 
4' 


1     sin  $  tan  0    cos  $  tan  0 
0  cos  $  —  sin  $ 

0     sin  $  sec  0     cos  $  sec  0 


(3.7) 


By  integrating  Equation  3.7,  the  time  history  of  the  Euler  angles  can  be  obtained. 

The  Euler  angle  method  has  one  drawback.  In  the  kinematic  differential 
equations  derived,  a  singularity  occurs  for  some  particular  value  of  either  $,  0,  or  ^. 
In  Equation  3.7,  this  singularity  occurs  at  0  =  ±7r/2,  which  means  the  that  the 
coordinate  transformation  is  useful  for  an  aircraft  in  horizontal  flight,  but  useless 
for  an  aircraft  which  requires  a  vertical  take-off  or  extended  periods  of  vertical  flight. 
Thus,  another  type  of  transformation  is  necessary  for  aircraft  that  spend  considerable 
time  flying  at  conditions  near  0  =  7r/2. 

The  first  alternative  is  simply  to  use  another  one  of  the  12  possible  Euler 
angle  transformations.  It  has  been  shown  that  the  rotation  matrix,  R,  is  made  up 
of  sequential  rotations  and  can  be  characterized  as  the  product  of  three  individual 
matrices  where 

/?-M,(^3)M/3(^2)Ma(^l),  (3.8) 

where  the  rotation  sequence  (a,  /?,  7)  represents  one  of  the  combinations  of  integers 

(1,2,3),(1,3,2),(1, 2,1), (1,3,1) 
(2,1, 3), (2,3,1), (2.1, 2), (2,3,2) 
(3,2,1), (3,1, 2), (3.2,3). (3,1, 3) 
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and  the  Mt{Oj)  are  the  rotation  matrices 


Miie)  = 


Ms{e)  = 


1         0  0 

0      COS  6      sin  6 
0    —  sin  ^    cos^ 

cos  6    0    —  sin  ^ 

0       1         0 
sin  6    0      cos  6 

cos  9      sin  6    0 

—  sin  6    cos  $    0 

0  0       1 


(3.9) 


One  can  see  that  by  substituting  in  (3,2,1),  and  $,0,  and  ^  for  ^i,^2i  and  ^3  re- 
spectively, the  matrix  shown  in  Equation  3.4  is  obtained.  Euler  angle  transformation 
matrices  for  each  of  these  combinations  have  been  calculated  and  are  tabulated  in 
[Ref.  Ka  83].  The  best  Euler  angle  rotation  sequence  for  an  aircraft  with  flight  at 
0  =  7r/2  was  selected  as  a  2-3-1,  or  Y-Z-X,  sequence.  This  sequence  will  allow  flight 
at  0  =  7r/2  with  no  corresponding  singularity  in  the  kinematic  differential  equations. 
Note  that  this  matrix  is  invertible  for  ^  7^  7r/2.  The  2-3-1  rotation  matrix,  R,  is 
described  in  terms  of  Euler  angles  as 

cos  0  cos  ^  sin  ^  —  sin  0  cos  ^ 

—  cos  0  sin  ^  cos  $  +  sin  $  sin  0      cos  ^  cos  0        sin  0  sin  ^  cos  0  +  sin  0  cos  0 
cos  0  sin  ^  sin  0  +  cos  $  sin  ^       —  cos^sin0     —  sin0sin  ^  sin  $  +  cos0cos  $ 

(3.10) 

The  kinematic  differential  equations  can  be  found  in  [Ref.  Ka  83]  as 

1 


0 


cos^ 


1     —  cos^sin^     sin  $  sin 'I' 
0  cos  $  -  sin  $ 

0      sin  $  cos  ^       cos  0  cos  ^ 


(3.11) 


These  equations  can  now  be  integrated  to  find  the  time  history  of  the  Euler  angles. 
2.     Quaternions 

Another  choice  for  the  expression  of  a  body's  spatial  orientation  is  the 
use  of  quaternions.  Quaternions  eliminate  the  disadvantage  of  the  singularity  in  the 
second  rotation  that  is  associated  with  the  Euler  angles.   Quaternions  have  been  in 
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use  for  quite  some  time,  having  been  discovered  by  Euler  in  a  searcli  for  complex 
numbers  [Ref.  Mo  84].  A  quaternion  is  defined  as  [Ref.  Mo  84]: 


q  =  l3o^il3,^jl32^ki3s, 


(3.12) 


where  the  parameters  are  represented  by  various  authors  as  S,  a,  b,  c  by  [Ref.  Ro  58], 

X,  <^,  r/,  and  (  by  [Ref.  Whi  59],  and  94,  ^i.  q2,  and  93  by  [Ref.  Sil  91]. 

The  components /?o,  ^1, /^2i  and /^3  are  real  numbers  and  the  terms  i,j,  and  k 

axe  defined  in  the  typical  manner  for  complex  numbers,  where 

P  =  —1  ij  =  —  j?  =  k 
p  =  -l  jk  =  -kj  =  i 
P  = -I     ki  =  -ik  =  j 

The  norm  of  a  quaternion,  q'q,  is  required  to  be  1: 


qq'  =  q'q  =  ^l^r  /3f  +  ^l  +  3l  =  1, 

since  q'  =  ^q  -  ?7?i  -  j^2-  kS^. 

It  can  be  shown  that  R  can  be  represented  as  follows: 


f,R  = 


ru 

r\2 

^13    ' 

^21 

r22 

^23 

''SI 

^32 

^33   . 

wh 


ere 


(3.13) 


(3.14) 


ru  =  2{M2^3o33) 

ri3  =  2{l5il33-i3oi32) 

r2i  =  2{i3,(32-l3o33) 

r22  =  3^  -  /??  +  /3|  -  3i  . 

r23  =  2(/?2/?3  +  ^0/?l) 

r3i  -  2(/9i^3  +  i5o^2) 

r32  =  2{M3-l3o3i) 

r33  =  /^o'-/^?-/^l+5| 
All  that  is  required  now  is  to  determine  the  kinematic  differential  equations  using  the 

quaternions. 
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By  expressing  subsequent  rotations  from  one  coordinate  system  to  another, 
where  /?'  orients  F^  to  F,  and  ^"  orients  F2  to  Fi,  an  algebraic  approach  [Ref.  Mo  84] 
can  be  used  to  relate  F2  to  F. 


R{/3)     =    R{l3")R{/3') 


(3.15) 


Now  the  ^'s  can  be  expressed  in  terms  of  /3"s  and  /?"'s  with  the  following  result 


fi  =  Rm^\ 


(3.16) 


Rm  = 


(3.17) 


where 

/^O      -A       -/^2      -/^3 
/9i  i5o         /^3      -/?2 

^2      -/?3  /3o  ^1 

By  regarding  the  second  rotation  in  Equation  3.1  Gas  infinitesimal,  the  following  result 
is  obtained 


^='^R{uj*)(3, 


where 


^  = 


'  k' 

r/?oi 

^ 
h 

,f3  = 

1^2 

k 

.03. 

and 


R{uj*)  = 
Equation  3.18  can  be  rewritten  as 


and  to*  = 


0  — u,'i  —uJ2  —<-03 

ix^l  U  (-J3  — UJ2 

OJ2  —'^3  0  Li>-[ 

U3  UJ2  —Wi  0 


1 


0 

p 

9 
r 


/?  =  -7?(/?)a*, 


(3.18) 


(3.19) 


(3.20) 


(3.21) 


and  in  this  form  can  be  integrated  to  give  the  time  history  of  the  orientation  of  the 
body. 
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Using  quaternions  has  the  following  advantages  over  Euler  angles  in  repre- 
senting spatial  orientation  of  a  rigid  body: 

•  Four  states  required  to  express  the  spatial  orientation. 

•  Requires  almost  30  %  fewer  calculations  [Ref.  Ro  58],  mainly  because  no  non- 
linear, trigonometric  equations  need  to  be  calculated. 

•  No  singularities  in  Equation  3.21  at  any  body  attitude. 

D.       DERIVATION  OF  EQUATIONS  OF  MOTION 

The  derivation  of  equations  of  motion  for  a  general  six  degree  of  freedom  airplane 
model  can  be  divided  into  two  parts.  The  first  part  is  simply  the  determination 
of  the  equations  of  motion  for  any  rigid  body  in  space.  It  is  dependent  only  on 
the  linear  and  angular  momenta  of  the  body.  The  second  part  is  the  calculation 
of  aerodynamic,  gravitational,  and  thrust  forces  on  the  airplane.  These  forces  are 
particular  to  a  certain  aircraft  and  in  general  can  be  represented  by  the  stability  and 
control  derivatives  described  later  in  the  thesis. 
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1.      Linear  Equations 

Equations  for  linear  motion  can  be  calculated  using  Newton's  law,  F  =  m  a. 
Since  most  of  the  sensor  information  available  for  feedback  to  the  control  and  nav- 
igation systems  is  available  in  the  [B]  coordinate  system,  the  terms  for  linear  ac- 
celerations, as  well  as  forces  and  moments,  will  be  expressed  in  the  body  coordinate 
system.  First  the  position  of  the  aircraft  e.g.  is  determined  as  ^' Pbo-  Then  Coriolis' 
theorem  is  applied  to  obtain  linear  velocities  for  the  aircraft.  Coriolis'  theorem  is 
then  reapplied  to  derive  the  expression  for  linear  accelerations.  Then 

""Vbo  =  "Pbo.  (3.22) 

Both  sides  of  Equation  3.22  are  premultiplied  by  fjR  to  get: 

yR  Vbo  =  irR  Pbo 

or 

^VBo  =  ^Pbo  (3.2.3) 

Now  consider  Coriolis'  theorem 

i  =  ^A  +  LjxA,  (3.24) 

at 

where  A  and  j-^A  use  the  notation  for  derivatives  previously  defined  in  Section  A.  The 

term  uj  x  A  represents  the  difference  between  the  relative  velocity  as  measured  from 

the  rotating  and  non-rotating  axes  [Ref.  Gre  88]. 

Equation  3.24  is  applied  to  ^vbo  in  Equation  3.23  to  get: 

^VBO  =  -T.^VBO  +  ^^B  X  ^VBO-  (3.25) 

at 
Newton's  law  can  now  be  written  as 

^F    =    m''a 

=    m  Vbo,  (3.26) 
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where  ^' F  is  the  total  external  force  applied  to  the  aircraft  e.g.    Equation  3.26  is 
premultiplied  by  yR  to  obtain  the  result: 

^F    =    mfjR^'hBo 

=    m^VBO.  (3.27) 


when  ^vbo  is  substituted  into  Equation  3.27,  the  final  result  for  ^F  is 

^F     =     m{-^VBO^^uJB  ^^VBO) 

=    m—^VBo-^rn^LjB^^VBO.  (3.28) 

at 

2.      Angular  Equations 

The  equations  for  angular  accelerations  are  derived  using  Euler's  law  for 
preservation  of  angular  momentum.  These  equations  are  also  derived  for  the  aircraft 
e.g.  by  applying  Coriolis'  theorem  to  the  equation  for  Euler's  law: 

''Lbo  =  ''Abo,  (3.29) 

where  ^' Lbo  is  the  angular  momentum  of  the  aircraft  and  ^^ Nbo  is  the  total  ex- 
ternal moment  applied  to  the  aircraft  e.g.  Euler's  law  can  be  rewritten  in  {B]  by 
premultiplying  Equation  3.29  by  fiR  to  get 

^Lbo  =  ?R"Nbo-  (3.30) 

Using  Coriolis'  theorem  in  Equation  3.24,  ^ Lbo  can  be  rewritten  as 

Lbo  =  -r   Lbo  +    ^B  X     Lbo-  (3.31) 

at 

The  angular  momentum,  ^ Lbo-,  is  defined  as  the  product  of  the  inertia  tensor,  Ib, 
and  the  body's  angular  velocity,  ^u,'s,  or 

L  —  iB      '^B  +  Ir      ^7?, 
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where  In  and  ^u^r  are  the  moment  of  inertia  and  the  angular  velocity  of  the  propeller, 
respectively.  When  this  term  is  substituted  into  Equation  3.31,  the  result  is 

d 


B 


LbO  =  J^Ub^'^B  +  iR^iOR)  +  ^UB  X  (/b^c^b  +  Ir^u;r), 


(3.32) 


where  Ib  is  the  inertia  tensor  for  the  aircraft  and  Ir  is  the  inertia  tensor  for  any 
significant  spinning  object  on  the  aircraft,  such  as  a  propeller,  turbine  ,  or  other 
rotating  disk.  The  term,  ^ojr,  is  the  angular  velocity  of  the  rotor,  expressed  in  {B}. 
We  can  carry  out  the  differentiation  in  Equation  3.32  to  get 


^Lbo  =  Ib-t^^b  +  Irt/u;r  +  ^ujB  X  {Ib^ujb  +  Ir^uJr). 
at  at 


(3.33) 


However,  since  d/dt{^iOB)  =  ^^b  and  d/dt{^u;R)  is  assumed  to  be  very  small.  Equa- 
tion 3.33  results  in 


Br 


B 


b 


Lbo  =  Ib  ^b  -\-    <^B  X  {Ib  ^b  +  Ir  ^r) 


Now  the  result  in  Equation  3.34  can  be  substituted  into  Equation  3.29: 


(3.34) 


B 


B 


B 


B 


Nbo  =  Ib  ^b  -\-    ^B  X  {Ib  u:b  +  Ir  ^r) 


(3.35) 


The  term  Ir^lor  can  be  disregarded  if  it  is  insignificant  compared  to  Ib  and     ljb 
[Ref.  Ros  79]. 

3.     State  Equations 

In  the  preceeding  sections,  kinematic  equations  for  the  motion  of  a  rigid 
body  were  derived  in  matrix  form.  These  equations  can  be  assembled  into  a  state 
space  representation  of  the  kinematic  equations  of  motion.  First,  Equations  3.28  and 
3.35  can  be  written  as 


Bf 


m 


i^VBo  +  rn  {^LCB  ^^VBo) 


Ib^^b  +  ^^B  X  {Ib^^b  +  Ir^^r)  . 


(3.36) 
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Equation  3.36  can  be  rearranged  to  yield 


Jt 


m  ^VBo 


-^LJB  X  {Ib^^b  +  Ir^^r)    +^A') 


(3.37) 


The  terms  on  the  left  hand  side  of  Equation  3.37  can  be  normalized  by  multiplying 
by  1/m  and  7g\  with  the  final  result  of 


di 


B 


VBO 


LOB 


-^uJB   X  ^VBO 


+  ^ 


(3.38) 


-Ib'^u:bx{Ib^u:b  +  Ir^^r)    +    Ib'^^' 
E.     EXTERNAL  FORCES  AND  MOMENTS 

Section  D.  gives  the  equations  for  kinematic  motion,  as  shown  in  Equation  3.38, 
for  any  rigid  body.  Now  it  is  necessary  to  distinguish  between  the  different  platforms 
to  be  modeled  in  order  to  give  an  accurate  representation  of  the  aircraft.  This  is 
achieved  by  computing  the  actual  ^F  and  ^A'  acting  on  the  aircraft.  These  forces 
and  moments  are  those  due  to  gravitational,  propulsive,  and  aerodynamic  effects, 

written  as 

F      _        Fqrav -\-    FpBop -^    Faero  /q  QQ^ 

B\j     -  B\r  I  fiv  (.o.-^yj 

1.      Aerodynamic  Forces  and  Moments 

The  aerodynamic  force  and  moment  terms  are  determined  by  using  a  first- 
order  Taylor  series  expansion  around  a  given  nominal  operating  point.  This  operating 
point  is  the  state  chosen  to  represent  the  aircraft's  flight  condition.  Each  term  in  the 
series  is  a  partial  derivative  of  ^F  or  ^jV  with  respect  to  the  aerodynamic  variables, 
ulU,  Q,  (3,  p,  g,  r  [Ref.  Sch  92,  Th  89]: 


Faero  =  SFj.>x'  -f  SF^'X  4-  SF^A  +  Fq. 


(3.40) 


Similarly,  moment  terms  can  be  written  as 


Naero  =  SN^'x'  +  SNr'X  +  SN^A  +  Aq, 


(3.41) 
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where  x'  is  given  by 


and  i'  is  given  as 


'  -\—   fi        pb_    qc_    rb_. 


i'  =  ['fi.a]. 


(3.42) 


(3.43) 


Control  inputs  are  represented  by  the  vector  A: 


A=   [8,,8r.6a\ 


(3.44) 


where  6e,  ^r,  and  8a  are  the  elevator,  rudder,  and  aileron  inputs,  respectively. 
Equations  3.40-3.44  can  now  be  combined  as  follows: 


w  p 


w^ 


=  qS{—/  +  —/  +  --A  +  Cfo). 
ox'         ox'         oA 


(3.45) 


w 


here  q   =    l/2pV^,   S   =   diag{5, 5, 5, 56, 5c,  56},  and   C  is  the  matrix  of  non- 
dimensional  stability  derivatives  differentiated  with  respect  to  the  terms  defined  in 


Equation  3.42,  3.43,  or  3.44.  |g  is  defined  as: 


Clv     Cl 


Cd 


Clc  Clj,  Cl^  Cir 

Cy0  Cy^  Cy^  Cy^  C)v 

Cd0  Cdo.  Cpp  Coq  Cor 

Ciu     Ci^  Ci^  Cip  Ci^  Ci^ 

r       c  c  c  c  c 

^mu       ^^vn.0  ^  ma  ^'ra-p  ^  niq  ^nir 

r       r  r  c  c  r 

^ni!        ^'na  ^Wq  ^"rip  '^'Uq  '^'tit 


jp  is  very  similar  to  ^,  except  that  only  the  a  and  /i  terms  are  normally  computed, 
leaving  a  6  x  2  matrix  rather  than  a  square  matrix.  The  term  |£  is  defined  as 


C 


Cr^        Cl 


Cys,  Cy,,  Cy^^ 

Cds,  Coi^  Cof^ 

^h,  ^'hr  ^'6a 

c  c  c 

^mi^  ^TTii^  ^msa 


^n/,.        ^'nt. 


^6e 


c 

^  riA 


23 


Cfq  is  defined  to  be  the  vector  of  steady  state  coefficients: 


Cfo  = 


Cdo 
Cyo 
Clo 

ClQ 
.    CnO 


representing  conditions  in  trimmed,  balanced  flight.  This  definition  is  similar  to  the 
definition  used  by  Roskam  [Ref.  Ros  79].  In  other  references,  the  term  Cfo  can 
refer  to  the  nominal  value  of  the  coefficient  at  a  =  0.  However,  in  the  Taylor  series 
expansion  it  is  more  natural  to  use  the  first  definition  of  Cfq:  therefore,  it  will  be 
used  in  the  following  derivation  and  modeling.  The  stability  and  control  derivatives 
are  usually  computed  in  the  so-called  wind  axis  coordinate  system.  The  wind  axis 
coordinate  system,  {H),  is  defined  as  the  coordinate  system  that  results  when  the 
xb  axis  is  aligned  into  the  relative  wind.  This  axis  will  not  be  aligned  with  the  body 
coordinate  system  since  the  aircraft  flies  with  an  angle  of  attack,  q,  and  can  have  a 
sideslip  angle,  l3.  The  transformation  from  {W]  to  {B}  is  performed  in  the  same 
fashion  as  the  Euler  angle  transformations  mentioned  earlier.   The  rotation  matrix, 


^rR,  is  a  function  of  a  and  /3,  and  is  expressed  as 


^■R 


cos  a  cos  (3    —  cos  asm  /S    —  sin  a 

sin  /?  cos  0  0 

sin  Q  cos /5     —  sinQsin/3      cosq 


(3.46) 


The  rotation  from  {W}  to  {B}  is  made  by  premultiplying  the  force  or  moment  vector 
by  ^/?.  Additionally,  since  the  lift  and  drag  are  defined  as  positive  along  the  negative 
zb  and  xb  axes,  we  define  Faero  and  Naero  as 


(3.47) 


'  -D  ' 

'   I 

\ero  = 

y 

and  Naero  = 

m 

-L 

n 

24 


In  order  to  write  Equation  3.45  in  state  space  form,  state  variables  must 
be  defined.  The  most  commonly  used  notation  to  use  for  the  state  vector  is  to  use 


X  = 


w 
V 

r 


(3.48) 


However,  the  terms  x'  and  x'  in  Equations  3.40  and  3.41  cannot  be  used  directly  as 
states.  Define 


M'  :  T  ->  x' 
M'  :i^  X 


(3.49) 


wnere 


and 


M'  =  diag{l/VT,l/V^T,l/VV,V2Vr,c/2VT,6/2VT} 


M'  = 


0    c/{2Vt)         0         0    0    0 
0         0         6/(2Vr)    0    0    0 


are  matrices  of  appropriate  dimensions.    The  complete  expression  for  ^ Faero  and 
^ ^AERO  can  now  be  written  as 


B 


S7V 


Faero 

AERO 


=  qS 


B 

W 


R      0 

B 
W 


0      S.i? 


{CFo  +  ^M'x  +  lgyi/'i  +  ^A}         (3.50) 


dx'  dx'  dA 

and  can  be  substituted  into  Equation  3.38. 
2.     Other  Forces  and  Moments 

In  addition  to  the  forces  and  moments  due  the  aerodynamics  of  the  aircraft, 
forces  and  moments  due  to  the  propulsive  and  gravitational  forces  must  be  considered. 
Gravitational  forces  acting  on  the  aircraft,  ^Fgrav,  can  be  found  by  premultiplying 
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Fgrav  by  the  appropriate  rotation  matrix,  where 

Fgrav  = 
Then 


0 

0 

mg 


Fgrav  =  vR  Fgrav-  (3.51) 

Propulsive  forces  and  moments,  ^ Frrop  and  ^ Nrrop,  are  computed  directly  in  [B] 


and  can  be  expressed  as: 


'F 


PROP  = 


and 


B 


NpRop 


Tx 

Ty 
Tz 


T, 

Tm 
Tn 


(3.52) 


(3.53) 


where  T,'s  represent  the  forces  or  moments  due  to  thrust.  Computation  of  propulsive 
forces  and  moments  depends  on  each  particular  engine  installation,  and  must  be 
determined  for  the  individual  aircraft  modeled. 

Equations  3.51,  2.,  and  2.  can  now  be  substituted  into  Equation  3.38: 


d 
di 


VBO 


B 


UJB 


B 


u^'fiX  0 

0  -^/fi^^c-fi  X  {^Ib^u^^b  +  Ir^u^^r)) 

-L        0 

TO 

0     «/b' 


B 


vbo 

!     , 

"  ^F 


+ 


(3.54) 


where 


Bf 


^K 


GRAV 
0 


+ 


^PROP 

B  ,\' 
l^PROP 


+ 


gi?      0 

0     g./? 


q'S  {  Cpo  +  fgM'x  +  IfA/'i  +  fgA  }  I  I .  (3.55) 


In  order  to  write  Equation  3.38  in  state  space  form,  the  terms  associated  with  x'  must 
be  collected  and  moved  to  the  left  hand  side,  along  with  the  other  time  deri\'ative 
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terms,  ^vbo  and  ^ujb-  Let 


w 


T  = 


0  ^rR 


and  Mj  = 


m       0 
0     ^1 


B 


then  the  complete  non-linear  equations  of  motion  for  any  aircraft  can  expressed  in 
state  space  form  as  follows  [Ref.  Th  89]: 


Jt 


B 


VBO 


B 


LOB 


=  X 


-1 


-^ujbX 
0 


0 

B  T    B. 


B 


-^/5H^u;5x(«/B«c^s  +  //?^c^fl)) 


MT'^TqS'-§^M' 


VBO 
^^B 


+  M/-' 


t^GRAX 

0 


^Fp 


syv 


ROP 
PROP 


St  +  ^^TqSiCf,  +  ^A) 


+ 
(3.56) 


and 


where 


Pbo  =  r-R   ^bo-, 


B, 


A  =  Sik^UB 


X  =  h-Mr^^TqS^M' 

ox' 


(3.57) 


(3.58) 


(3.59) 


P  is  the  position  vector  of  the  aircraft,  and  A  is  the  matrix  of  kinematic  differential 
equations  based  on  Euler  angles  or  quaternions. 
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IV.  COMPUTER  MODELING  OF  THE 
AIRCRAFT  EQUATIONS  OF  MOTION 

A.     IMPLEMENTATION  OF  EQUATIONS 
1.      Procedure 

The  ARGD/Archytas  and  the  Bluebird  aircraft  models  require  different 
non-linear  equations  of  motion.  This  difference  is  due  to  the  unique  nature  of  each 
aircraft.  In  the  case  of  the  AROD,  the  angular  momentum  of  the  propeller  is  a  major 
factor  to  be  considered,  while  the  aerodynamics  effects  in  a  hover  are  negligible. 
The  Bluebird  on  the  other  hand  is  a  conventional  aircraft,  and  exhibits  the  opposite 
characteristics.  Angular  momentum  from  the  propeller  is  small  and,  the  aircraft 
requires  the  stability  and  control  derivatives  associated  with  aerodynamic  flight.  All 
the  equations  were  implemented  in  a  systematic  manner  following  the  same  general 
approach  for  either  aircraft  model.  The  commercial  products  MatLAB  and  SiMULINK, 
©  1990-1992  the  Mathworks,  were  chosen  for  the  modeling\  mainly  due  to  the  ease 
of  expressing  matrix  equations.  The  model  was  constructed  using  the  following  steps. 

•  Kinematic  equations  of  motion  were  coded. 

•  Gravitational  forces,  with  the  direction  cosine  matrix  represented  by  Euler  an- 
gles were  added  to  the  model 

•  Stability  and  control  derivatives  were  included  in  the  model,  as  well  as  engine 
thrust,  as  appropriate  for  the  aircraft  being  modeled. 


'All  code  is  listed  in  APPENDIX  C. 

28 


•  Engine  and  actuators  were  added  to  the  model. 

•  Sensors  for  control  systems  were  added  to  the  model. 

In  the  first  stage  of  modeling,  analytic  linearization  was  also  carried  out  to 
verify  the  computer  calculations.  This  was  done  by  analytically  linearizing  the  matrix 
formed  from  the  six  non-linear  equations  governing  the  kinematic  motion  of  the  body. 
Nominal  values  were  substituted  in,  and  the  results  compared  to  the  trimmed  and 
linearized  values  obtained  from  the  SiMULINK  program^.  The  next  step  required  the 
addition  of  gravitational  terms,  and  analytic  linearization  was  still  manageable.  The 
linearized  matrix  now  included  nine  equations  and  nine  states  with  the  Euler  angle 
direction  cosine  matrix,  DCM,  and  ten  equations  with  ten  states  for  the  quaternion 
DCM.  Nominal  values  for  each  case  were  again  substituted  into  the  linearized  ma- 
trix and  compared  to  the  linearized  model  derived  from  SiMULINK.  The  inclusion  of 
the  stability  and  control  derivatives  and  the  thrust  terms  presented  a  problem  that 
was  much  too  cumbersome  to  linearize  analytically.  Verification  of  the  data  at  this 
stage  was  accomplished  by  direct  comparison  of  the  dimensional  derivatives  result- 
ing from  numerical  linearization  of  the  plant.  In  the  case  of  the  AROD  the  results 
were  compared  with  data  published  by  Sandia  Labs  [Ref.  Wh  91];  for  the  Bluebird, 
eigenvalues  were  computed  and  then  compared  to  eigenvalues  obtained  by  classical 
analytic  methods  [Ref.  Sch  92]  and  [Ref.  Ros  79]. 
2.      AROD  Equations 

The  AROD  differs  from  the  Bluebird  primarily  in  that  the  aerodynamic 
forces  and  moments  are  negligible  while  the  craft  is  hovering.  The  powered  lift  does 
present  some  special  problems,  the  predominant  difficulty  being  the  gyroscopic  cou- 
pling of  the  spinning  propeller.  Another  important  consideration  is  the  moment  due 


^Data  is  tabulated  in  APPENDIX  B. 
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to  the  swirl  effect  of  the  air  from  the  propeller  striking  the  vanes  mounted  aft  of  the 

propeller.  These  forces  and  moments  are  computed,  then  substituted  in  for  ^ F  and 

^N  in  Equation  3.38. 

a.      Control  Forces  and  Moments 

All  of  the  applied  forces  and  moments  in  the  AROD  are  due  to  the 

four  control  vanes  mounted  aft  of  the  propeller.  The  vanes  can  be  moved  in  different 

combinations,  as  discussed  in  later  in  this  chapter,  in  order  to  maneuver  the  AROD. 

The  forces  and  moments,  ^ F  and  ^ N ^  acting  on  the  AROD  are  computed  from  a 

Taylor  series  expansion  around  a  nominal  hover  point.  Since  aerodynamic  terms  are 

negligible, 

r  Bp^ 1  ^c 

(4.1) 


^iV, 


Control 
Control 


= '"■^«£^' 


where 


A  -  [<5e,  6,,  S,]' 


•  9.  =  i/2/>v; 


/2 


•  Vi  is  the  induced  velocity  through  the  propeller  [Ref.  Pro  90]  and 
V?  =  T/{2Ap). 

•  A  is  the  inlet  area,  where  A  =  3.14  p. 

•  R  is  the  propeller  radius,  where  R  =  l.O  p. 

Notice  no  aerodynamic  forces  act  on  AROD  due  to  the  movement  of  the  elevator,  rud- 
der, or  aileron  controls.  Again,  this  is  because  the  model  of  interest  is  only  designed 
to  fly  in  a  stable  hover.  The  other  forces  and  moments  involved  in  these  calculations 
are  due  to  gravitational  and  propulsive  action,  and  are  described  next. 
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b.      Additional  Forces  and  Moments 

The  gravitational  force  term,  ^Fgrav,  is  computed  in  the  {B}  coor- 
dinate system  with  the  rotation  matrix,  qR.  For  the  AROD,  gi?  is  determined  by 
the  2-3-1  Euler  angle  rotation  sequence  defined  in  Equation  3.10  and  is  written  as 


Fgrav  =  uR    Fgrav  where     Fgrav  = 


0 

0 

mg 


(4.2) 


where 


GRAY     =    rng 


—  sin  0  cos  ^ 

sin0  sin  ^  cos$  +  sin$  cos  0 

—  sin  0  sin  ^  sin  $  +  cos  0  cos  $ 


(4.3) 


The  forces  and  moments  due  to  the  propeller  thrust  were  determined 
experimentally  [Ref.  Sto  93]  and  are  discussed  in  a  later  section.  For  the  AROD,  the 
propulsive  force,  ^ Frrop,  is  acting  completely  along  the  xb  axis 


B 


FpROP  — 


Ftx 
0 
0 


(4.4) 


where  Fj^  is  the  total  thrust,  determined  as  a  function  of  rpm. 

The  moment  resulting  from  thrust  is  ^ Nrrop  and  is  given  by  the 
vector  equation 


B 


NpRop  = 


h 
0 
0 


(4.5) 


where  It  is  the  rolling  moment  due  to  thrust.  This  rolling  moment  is  due  to  the 
swirl  of  the  air  as  it  leaves  the  propeller  and  strikes  the  control  vanes.  This  term  was 
determined  experimentally  as  a  function  of  thrust,  and  is  also  discussed  further  in  a 
subsequent  section. 

When  all  the  terms  are  collected,  the  total  force  and  moment  applied 
to  AROD  can  be  expressed  as 


Bp 


=  {   qAR^A^ 


syv 


FpROP 
PROP 


+ 


B 


Fgrav 
0 


(4.6) 
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The  force  and  moment  terms  in  Equation  4.6  can  be  substituted  into 
Equation  3.38  with  the  resulting  state  space  equation  given  by 


d 
It 


B 


VBO 


■     uu'B    X       VbO 


■1   B 


-Iq'  ""lob  X  {Ib"oJb  +  Ir"uJr) 


+ 


//m      0 
0       Ib' 


qAR^A^ 


B 


^Nf 


+ 


t^GRAV 
0 


(4.7) 


FpROP 

hnop 
Equation  4.7  can  now  be  written  in  the  state  space  form  and  programmed  using 

SiMULINK 


d 

'  ^VBO  ' 

f 

"    -^C^B 

dt 

^UJB 

n 

0 

■  //m      0 

0       Jb    \ 

0 

B  T    B, 


-^lB\''uJBxriB''u;B  +  lR''u,>n) 
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VBO 
lJb 
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^N 


FpRop 

PROP 


+ 


Fgrav 
0 


•  (4.8) 


3.     Model  Validation 

a.      Kinematic  Equations 

As  discussed  in  the  beginning  of  the  chapter,  the  first  step  in  the 
AROD  model  validation  was  to  linearize  the  non-linear  equations 


dt 


^vbo    =     l/m{-^UBX  ^vbo) 


—    tUB     -       7^  (-   u;bX(    1r  u;r+    1b   ^b), 


(4.9) 
(4.10) 


with  the  result  given  by^: 
d 


dt 


Sv 


—IjJq  X  VqX 

0        -^Ib\^Ir^^r  X  -(-'ox)^/b  +  ^/b-'ox) 

0 


+ 


-u7o  X     Ir 


Sv 


S^'R.   (4.11) 


The  nominal  values  for  the  AROD  hover  operating  point, 

10 

0 

0 

0 

0 

0 


X  = 


(4.12) 


''See  APPENDIX  A  for  a  description  of  the  linearization  process 
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can  be  substituted  into  Equation  4.11  to  get 


r  0  0  0 

0 

0 

0       1 

0    0    0 

0 

0 

10 

d 

Sv 

0    0    0 

-10 

0 

0 

Sv 

dt 

Su! 

0    0    0 

0 

0 

0 

Suo 

0    0    0 

0 

0 

-1.615 

0    0    0 

0 

1.606 

0 

(4.13) 


where  the  linearized  matrix  associated  with  S^ji  is  zero,  since  —loq  x  ^ If{  =  0.  These 
values  match  the  values  obtained  by  linearizing  the  model  using  SiMULINK. 

The  inertial  cross  coupling  from  the  propeller  is  evident  from  the  lin- 
earized results.  The  values  of  —1.6152  sec~^  and  1.6062  sec~^  are  defined  as  pitching 
moment  due  to  yaw  rate,  rrir,  and  yawing  moment  due  to  pitch  rate,  n,,  respectively. 
The  gyroscopic  coupling  is  demonstrated  by  putting  a  step  moment  input  of  1  second 
duration  along  the  y  axis  and  observing  the  time  history  of  the  angular  rates  q  and 
r  as  shown  in  Figure  4.1.  The  AROD  has  a  tendency  to  spin  in  the  axis  orthogonal 
to  the  input  torque  as  shown  by  the  motion  r  about  the  z  axis. 
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Figure  4.1:  Gyroscopic  Motion  of  AROD 
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b.      Gravitational  Forces 

The  next  step  in  modeling  the  AROD  was  to  include  gravitational 
effects  into  the  model.  This  required  the  expression  for  ^Fgrav,  given  by  Equa- 
tion 4.3,  to  be  coded  into  a  Matlab  function.  This  function  block  was  then  added 
to  the  SiMULINK  model,  as  shown  in  Figure  4.2.  The  equation  to  be  analyzed  at  this 
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Figure  4.2:  Modeling  of  Gravitational  Effects  for  AROD 
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where  5(A)  represents  the  kinematic  differential  equations,  defined  in  Equation  3.11. 
Notice  that  since  gravity  acts  through  the  e.g.,  no  ^Ngrav  term  exists.  Next,  Equa- 
tion 4.14  was  linearized  both  analytically  and  with  SiMULINK.  The  analytic  result, 
given  in  state  space  form,  was 
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Slor. 


(4.15) 


In  Equation  4.15,  the  functions  /(A),G(A),  and  /?(A,u;)  are  partial  derivatives  of 
Equations  4.3  and  3.11  with  respect  to  A,  or  d/dA.  The  function  /(A)  is  the  linearized 
expression  for  ^ Fqrav  given  by  Equation  4.3,  where 


/(A)- 


d  ^E 


dA 


GRAY   _      _d_ 

m       ~^  dA 


—  sin  0  cos  ^ 

sin  0  sin  ^  cos  $  +  sin  $  cos  0 

—  sin  0  sin  ^  sin  $  +  cos  0  cos  ^ 


(4.16) 


or 


0  —  cos0cos^  sin0sin^ 

—  sin0sin^  sin$  +  cos$  COS0        cos0  sin'P  cos<^  —  sin<I>  sin0       sin0cos^cos$ 

—  sin$cos0  —  sin0sin^  cos$    —  cos*l>  sin0  —  cos0sin'I' sinO    —  sin0cos^  sin$ 

(4.17) 

where  /(A)  is  evaluated  at  the  nominal  condition,  Aq.  The  function 


B 


G{A)  =  d/duj{S{AfiOB 


is  derived  by  linearizing  the  kinematic  differential  equations  in  Equation  3.11.  The 
function  h{A,uj)  =  d/dA{S{A)^UB)  is  not  presented  since  for  u-'o  =  0,  as  in  steady 
state  cruise,  /i(A,u;)  =  0.  The  matrix  6'(A)  is  derived  from  the  following  equation 
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(4.18) 


(4.19) 
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These  matrices  were  assembled  as  in  Equation  4.15  and  were  evaluated  for  the  nominal 
flight  condition.  The  values  for  xq  for  this  flight  condition  were  given  by 

10 
0 
0 
0 
Xo=       0       .  (4.20) 

0 
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0 
After  substitution  into  Equation  4.15,  the  result  for  the  linearized  equations  of  motion 

was 
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(4.21) 


This  result  was  essentially  indentical  to  the  results  obtained  by  using  the  trim  and 
linearize  functions  in  SiMULINK.  The  gravitational  effects  acting  on  the  aircraft  were 
examined  by  running  a  simulation  of  the  non-linear  model  and  comparing  the  results 
to  those  from  Equation  4.11.  The  expectation  is  that  a  body  falling  in  earth's  gravity 
will  experience  an  acceleration  of  32.174//5'^  and  as  shown  in  Figure  4.3,  this  expec- 
tation is  realized  in  the  model.  In  Figure  4.3,  at  the  end  of  the  10  second  simulation, 
the  vertical  velocity,  u  in  this  case,  is  ~  320//s,  as  was  predicted. 
c.      Additional  Forces  and  Moments 

The  last  step  in  modeling  the  AROD  was  to  include  the  forces  and 
moments  due  to  propulsive  and  control  action  that  act  on  the  aircraft.  The  data  used 
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Figure  4.3:  Velocity  of  AROD,  Gravitational  Effects  Included 

in  modeling  these  control  forces  were  collected  experimentally  [Ref.  Sto  93]  and  then 
curve  fitted  to  a  function.  The  measurements  of  rolling  moment,  thrust,  rpm,  and 
vane  position  were  taken  for  various  configurations.  The  data  collected  were  then 
reduced  and  the  required  characteristics  were  computed.  The  accuracy  of  the  AROD 
model  depended  on  very  accurate  modeling  of  thrust  as  a  function  of  rpm,  moment  as 
a  function  of  thrust  produced,  and  moment  produced  by  deflecting  the  control  vanes 
in  the  different  combinations. 

Thrust  and  rolling  moments  were  measured  directly  at  different  power 
settings  ranging  from  3000  rpm  to  7600  rpm,  with  a  power  setting  of  ~  6400  rpm 
giving  a  thrust  of  ~  90/6/.  This  thrust  is  aprroximately  the  force  required  to  maintain 
a  hover  for  the  basic  AROD  configuration.  Figure  4.4  shows  the  linear  curve  fit 
through  the  thrust  vs.  rpm  data.  The  curve  was  fit  using  data  from  5000  rpm  to 
7600  rpm,  as  this  was  expected  to  approximate  the  normal  operation  range  in  flight. 
The  thrust  and  moment  data  are  plotted  in  Figure  4.5  along  with  the  line  fit  through 
those  data  points.  The  best  fit,  by  least  squares,  for  the  data  in  Figure  4.4  was  given 
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Figure  4.4:  Thrust  vs.  RPM  for  AROD 
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Figure  4.5:  Thrust  and  Moment  Data  for  AROD 


by 


Ft,  =  0.0291  Srpm  -104.7, 


(4.22) 


where  Srpm  represents  the  rpm  at  a  given  throttle  setting.  The  best  fit  for  the  data 
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in  Figure  4.5  was  given  by 

Nt  =  -0.0542  Ft,  -  0.9138.  (4.23) 

These  equations  were  then  used  in  modeling  the  engine  response  to  a  rpm  setting. 

The  engine  itself  was  modeled  using  a  second  order  transfer  function 
from  the  servo  position  to  Srpmi  based  on  the  Sandia  Labs  [Ref.  Wh  87]  model.  The 
transfer  function  for  the  engine  was  given  as 

^rpm  =      2^^^"^^        2^T  (4.24) 

where  Ke  =  900,  u^  =  5rad/sec,  (  =  1.0,  and  Sj  is  the  throttle  servo  position.  Since 
the  actual  throttle  position  is  set  via  a  radio  link  and  tests  have  not  been  set  up  to 
model  the  response,  it  was  ignored  in  the  model. 

The  engine  servo  could  be  modeled  and  a  transfer  function  from  com- 
manded input  to  servo  position  was  determined.  This  transfer  function  also  was 
determined  by  Sandia  Labs  in  the  original  AROD  work  [Ref.  Wh  87]  to  be 

^^  =    2^of^    j_    2^^  (4-25) 

where  u;„  =  2Qrad/sec^  (  =  0.6,  and  uj  represents  the  commanded  input  to  the  servo. 
It  should  be  noted  that  Equation  4.25  was  also  used  for  the  control  vane  servos,  since 
all  the  servos  in  the  AROD  were  identical. 

In  order  to  model  the  moments  due  to  control  surface  deflection, 
^Control:  computation  of  control  vane  effectiveness  was  necessary.  Again,  the  data 
gathered  by  [Ref.  Sto  93]  was  used  to  compute  vane  effectiveness  for  the  AROD 
model.  Figure  4.6  shows  the  moment  data  for  75%  thrust.  This  power  setting  cor- 
responds closely  to  the  thrust  required  to  maintain  hovering  flight.  A  line  was  fitted 
through  the  data  points,  and  the  slope  of  this  line  was  used  to  determine  the  rolling 
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Figure  4.6:  Moment  Due  to  Vane  Deflection 


moment,/,  for  a  given  vane  deflection.  The  equation  for  the  curve  was 


/  =  0.5523  V^v'G- 5- '248 


(4.26) 


for  6580  rpm,  where  Vavg  was  the  average  vane  deflection,  1/4  ^f_j  \],  in  degrees. 
This  averaging  was  required  since  the  individual  vanes  were  at  slightly  different  posi- 
tions with  respect  to  the  commanded  position.  Computing  the  vane  effectiveness  for 
roll  was  done  by  differentiating  with  respect  to  the  vane  deflection. 


/.    - 


dl 


0.5523 


fi  -  Ibf 


(4.27) 


OVavg  deg 

The  rolling  moment  should  be  non-dimensionalized  for  use  in  the  equations  of  motion 
given  by  Equation  4.8.  This  was  done  by  applying  the  following  equation 


Ci 


(4.28) 


l/2pV256 

To  use  the  measurable  quantities  available,  redefining  the  terms  in  Equation  4.28  was 
necessary.  The  representative  area,  S,  is  defined  as  the  inlet  area  to  the  duct,  k.  The 
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characteristic  length,  b,  is  defined  as  the  propeller  radius,  R.  The  velocity  term  was 
defined  as  induced  velocity,  V,.  With  these  terms  the  non-dimensional  coeflficient  C)^ 
could  be  calculated  as 

and  was  computed  to  be  C/^  =  1 .438 /rad.  The  rolling  moment  was  measured  during 
the  testing  and  was  easily  non-dimensionalized  in  a  form  suited  to  computer  modeling. 
However,  no  measured  data  existed  for  the  pitching  and  yawing  moments.  This 
required  the  estimation  of  pitch  and  yaw  vane  effectiveness  by  a  simple  ratio  technique. 
First,  it  was  assumed  that  the  vane  effectiveness,  /^^  for  two  vanes  (2V),  would  be 
exactly  twice  the  /^^  for  one  vane,  and  1/2  the  /^^  for  four  vanes,  (4V').  This  would 
make 

oVavg  2  oVavg 


or  numerically. 


(#-h.  =  0.2762:^  (4.31) 

oVavg  deg 


Now,  the  ratio  of  the  dimensional  derivatives  to  the  moment  arms  was  set  up  as 

/       dl       \  f      dm      \ 

^  dVAVG  ''^^'   ^   ^dVAVc'  U  32) 

where  l^  was  the  distance  from  the  e.g.  along  the  x-axis,  to  the  midspan  of  the  control 
vane,  a  distance  of  9.0  in.  The  distance  ly  was  the  distance  from  the  e.g.  along  the 
y-axis,  to  the  1/4  chord  of  the  control  vane,  a  distance  of  15.43  in.  The  calculation 
was  performed  and  the  result  for  the  pitching  moment  was 

^^     -0.4735=^^  (4.33) 


dVAVG  deg 

This  was  non-dimensionalized  using  Equation  4.28  with  the  result  of  C^^    =  1 .2332rac?~^ 
Moreover  because  of  the  symmetrical  design  of  AROD,  Cmt^  =  C'n<^-  The  results  for 
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vane  effectiveness  are  important  to  a  high  fidelity  model  and  are  presented  along 
with  other  relevant  data  in  Table  4.1.  The  non-dimensional  derivatives  in  Table  4.1 

TABLE  4.1:  NON-DIMENSIONAL  DERIVATIVE  DATA  FOR  AROD 


rad  ^ 

deg-' 

Positive  Vane  Deflection 

(^'.. 

1.438 

0.0251 

V'l  -f  V2  +  V3  +  V4 

^'rn6^ 

-1.233 

-0.0215 

V4  -  V2 

^rif^ 

-1.233 

-0.0215 

V3  —  Vi 

were  substituted  into  the  term  [dCp/dA)  in  Equation  4.7.  Written  as  a  matrix,  the 
non-dimensional  derivatives  are 


dCp 
dA 


A  = 


1 

mi. 

0 

0 

'  s. 

0 

Cni^ 

0 

Sr 

0 

0 

Cu. 

Sa 

(4.34) 


Equations  4.22,  4.23,  and  4.34  were  added  as  a  separate  block  in  the  model,  resulting 
in  a  block  diagram  shown  in  Figure  4.7. 
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Figure  4.7:  Non-Linear  Equations  of  Motion  Model 
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In  the  SiMULINK  linearization  results,  no  change  was  expected  in  the 
A  matrix  since  no  aerodynamic  terms  were  added  to  the  model.  Rolling  motion  is 
shown  in  Figure  4.8  as  an  example  of  the  unstable  nature  of  the  AROD  without  a 
stability  augmentation  system. 
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Figure  4.8:  Rolling  Motion  For  Complete  Non-Linear  Equations 


4.      Bluebird  Equations  of  Motion 
a.      Kinematic  Equations 

Modeling  the  kinematic  equations  of  motion  for  Bluebird  was  accom- 
plished using  the  same  procedure  as  was  used  in  the  AROD  modeling.  A  slight  dif- 
ference arose  in  the  derivation  of  the  angular  acceleration  equation  for  Bluebird  since 
the  term,  Irur,  for  the  angular  momentum  due  to  propeller  rotation  is  negligible. 
Thus,  the  following  equations  were  used  in  the  first  stage  of  modeling. 

d_ 
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B. 


B  J    B 


dt 


^S      =      "IE   i-^B  ^"Ib"^b). 


(4.35) 
(4.36) 
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These  equations  were  coded  into  a  MaTLAB  function  and  placed  into  a  block  diagram, 
shown  in  Figure  4.9.   Equations  4.35  and  4.36  were  linearized  analytically'*  with  the 
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Kincmatjc 


Figure  4.9:  Block  Diagram  of  Kinematic  Equations  of  Motion 


resulting  state  space  equation 


It 


Sv 


0 


bv 

buJ 


(4.37) 


This  equation  was  evaluated  at  the  nominal  flight  conditions,  determined  by  the 
typical  cruise  condition  for  the  Bluebird  aircraft.  A  state  vector  for  the  nominal 
flight  conditions  is  given  by 


2-0 


72.9954  fis 

0 
6.6757  fIs 

0 

0 

0 


(4.38) 


''See  APPENDIX  A  for  the  details 
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These  values  of  xq  were  substituted  into  Equation  4.37: 


d_ 
di 


Slj 


0  0  0       0  -6.675          0 

0  0  0  6.675  0  -72.995 

0  0  0       0  72.995           0 

0  0  0       0  0               0 

0  0  0       0  0               0 

0  0  0       0  0               0 


6v 


(4.39) 


These  results  were  in  complete  agreement  with  the  data  obtained  by  trimming  and 
linearizing  the  non-linear  model  with  the  SiMULINK  program.  The  comparison  be- 
tween Equation  4.39  and  Equation  4.13  shows  the  absence  of  any  angular  rate  cross 
coupling  in  Bluebird.  The  absence  of  the  cross  coupling  terms  results  from  the  choice 
to  ignore  the  angular  momentum  contribution  from  the  propeller. 
b.      Gravitational  Forces 

After  the  basic  kinematic  equations  Equation  4.35  and  4.36  were  put 
into  block  diagram  form,  it  was  an  easy  matter  to  include  additional  blocks.  The 
next  block  was  a  function  block  including  the  ^ Fgrav  terms.  The  model  with  the 
gravitational  terms  included  is  shown  in  Figure  4.10.  The  equations  of  motion  to  be 
modeled  at  this  stage  were 
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^A    =    S[\fu:B 


(4.40) 
(4.41) 
(4.42) 


where  5(A)  is  the  set  of  kinematic  differential  equations  based  on  the  3-2-1  Euler 
angle  rotation  in  Equation  3.7.  These  equations  were  then  linearized  analytically, 
with  the  result 
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Figure  4.10:  Gravitational  Forces  Model 

where  the  terms  /(A),  G(A),  and  /i(A,a;)  are  based  on  the  3-2-1  Euler  rotation  se- 
quence. The  linearization  is  done  in  the  same  manner  as  was  shown  in  Equations  4.16, 
4.17,  and  4.19  resulting  in 


and 


Note  that 
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MA,  c.-)-^5(A)^u;b|o  =  0, 
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(4.45) 
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since  uq  =  0.  When  Equation  4.43  is  evaluated  at  the  nominal  condition  xq  given  by 

72.9954  f/s 

0 
6.6757  f/s 

0 
xo=  0 

0 

0 
0.0912  rad 

0 


the  resulting  equation  is 
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(4.46) 


In  this  instance,  as  well,  the  results  of  the  analytic  linearization  in  Equation  4.46 
agree  very  closely  with  the  computed  results. 

A  non-linear  simulation  of  the  system  in  Figure  4.10  should  show  an 
increasing  downward  velocity  due  to  the  gravitational  effects  of  ^ Fgrav-  One  would 
also  expect  to  see  decreasing  forward  velocity  due  to  the  "drag-like"  term  that  arises 
with  the  introduction  of  the  angle,  0.  This  plot  is  shotvn  in  Figure  4.11. 
c.      Aerodynamic  Forces  and  Moments 

Completion  of  the  Bluebird  equations  of  motion  model  required  the 
modeling  of  the  aerodynamic  forces  and  moments  acting  on  the  aircraft.  A  simple 
engine  model  was  also  developed.  No  analytic  linearization  was  performed  at  this 
stage  due  to  the  increased  complexity  of  the  model.  Verification  of  the  computer 
results  was  accomplished  by  comparing  the  modes  and  eigenvalues  of  the  computed 
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Figure  4.11:  Gravitational  Effects  on  Velocity 

plant  with  those  resulting  from  substituting  the  stability  and  control  derivatives  into 
equations  developed  in  [Ref.  Sch  92]. 

The  aerodynamic  forces  and  moments  as  described  in  Equation  3.50 
were  coded  as  a  MaTLAB  function,  then  included  as  a  block  in  the  model  shown  in 
Figure  4.12. 

Next,  it  was  necessary  to  premultiply  all  the  blocks  by  \~\  since 
this  added  the  effects  of  the  d  derivatives.  Now  the  important  task  was  to  calculate 
the  stability  and  control  derivatives  using  the  general  reference  for  the  estimation  of 
non-dimensional  derivatives,  DATCOM  [Ref.  USAF  60].  The  stability  and  control 
derivatives  were  computed  based  the  aircraft  geometry  and  the  control  surfaces.  These 
values  are  tabulated  in  Table  4.2  and  in  Table  4.3  where  the  non-dimensional  force  is 
listed  on  the  left  side,  and  the  particular  derivative  is  determined  using  the  top  row. 
For  example,  Coa  —  0.188  using  Table  4.2.  The  terms  in  Cfq  were  also  estimated  to 
be  Clo  =  0.385  and  Cdq  =  0.03,  with  all  other  terms  equal  to  zero  since  the  aircraft 
is  in  straight  and  level  flight  at  this  trim  condition. 
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Figure  4.12:  Full  Non-Linear  Equations  of  Motion  Model 


TABLE  4.2:  NON-DIMENSIONAL  STABILITY  DERIVATIVES 


u 

(i 

Q 

P 

q 

r 

Q 

Cd 

0 

0 

0.188 

0 

0 

0 

0 

Cy 

0 

-0.31 

0 

0 

0 

0.0973 

0 

Cl 

0 

0 

4.22 

0 

3.94 

0 

1.32 

Ci 

0 

-0.0597 

0 

-0.363 

0 

0.100 

0 

Cm 

0 

0 

-1.163 

0 

-11.77 

0 

-4.70 

Cn 

0 

0.0487 

0 

-0.0481 

0 

-0.0452 

0 

The    FpROP  was  based  on  estimating  engine  thrust  for  a  4  HP  engine 
and  a  propeller  efficency,  rjp,  of  0.65.  Thrust,  To  was  estimated  using  the  equation 

bbOrjpHPpn 


T  = 


(4.47) 


^^0         Po 

where  /)„  is  the  density  at  the  operating  altitude,  Pq  is  the  sea  level  density,  and  f^o  is 
the  velocity  of  the  aircraft  in  f/s.  The  result  was  an  engine  thrust  of  Tq  =  19.5/6/,  for 
a  density  ratio  of  1.  This  could  be  directly  factored  into  the  equations  of  motion  as 
TqSj,  where  Sj  w'as  arbitrarily  scaled  from  zero  to  one.  Sj  =  0  represents  zero  thrust 


49 


TABLE  4.3:  NON-DIMENSIONAL  CONTROL  DERIVATIVES 


Se 

Sr 

Sa 

Cd 

0.065 

0 

0 

Cy 

0 

0.0697 

0 

Cl 

0.472 

0 

0 

Ci 

0 

0.0028 

0.265 

Cm 

-1.410 

0 

0 

Cn 

0 

-0.0329 

-0.0347 

and  St  =  I  represents  maximum  thrust. 

The  preceeding  values  were  substituted  into  the  appropriate  Matlab 
functions  and  the  entire  equations  of  motion  model  was  trimmed,  linearized,  and 
then  compared  with  the  analytic  results  based  on  classical  techniques  for  determining 
eigenvalues  and  eigenvectors.  The  eigenvalues  are  compared  in  Table  4.4.  The  results 
from  the  computed  eigenvalues  is  very  close  to  the  eigenvalues  derived  by  analytic 
methods. 

TABLE  4.4:  COMPARISON  OF  EIGENVALUES  FOR  BLUEBIRD 


MODE 

COMPUTED 

ANALYTIC 

LONGITUDINAL 

Phugoid 

-0.0191  ±0.4963j 

-0.0473  ±  0.4940j 

Short  Period 

-3.9833  ±3.5521; 

-4.0034  ±  3.5462; 

LATERAL 

Dutch  Roll 

-0.5285  ±  3.6346; 

-0.522  ±3.6194; 

Short  Period 

-5.5629 

-5.6654 

Spiral 

+0.0420 

+0.0420 

B.     VALIDATION  OF  AN  INDEPENDENT  CASE 

Although  verification  of  the  model  was  accomplished  at  each  stage,  comparison 
of  the  results  from  the  numerical  linearization  with  linearized  results  from  an  inde- 
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pendent  source  was  still  necessary  before  the  model  could  be  considered  completely 
reliable.  The  test  case  selected  was  the  Cessna  172  documented  in  [Ref.  Ros  79]  as 
airplane  A.  The  tabulated  non-dimensional  derivatives  and  given  flight  conditions 
were  used  as  inputs  to  the  non-linear  model.  The  model  was  then  trimmed  at  the 
specified  flight  condition  of  Vj  =  219  f/s  and  0o  =  0.  Using  the  resulting  state  and 
input  vector,  the  model  was  linearized  around  thesed  nominal  conditions.  With  the 
linearized  plant,  eigenvalues  were  determined  and  compared  to  those  tabulated  for 
airplane  A  (see  Table  4.5).  Very  little  difference  between  the  eigenvalues  is  seen  in 
the  longitudinal  modes.  Slightly  greater  differences  are  noticed  when  comparing  the 
lateral  modes,  but  these  differences  are  still  fairly  small.  Another  very  good  method 

TABLE  4.5:  EIGENVALUE  COMPARISON  FOR  CESSNA  172  TEST 
CASE 


Mode 

Numerical 

Independent 

Longitudinal 

Short  Period 

-4.130  ±4.3895j 

-4.130  ±4.390j 

Phugoid 

-0.0209  ±0.1794; 

-0.02092  ±0.1 797; 

Lateral- Directional 

Dutch  Roll 

-0.6947  ±  3.3080; 

-0.6858  ±  3.306; 

Spiral 

-12.4309 

-12.43 

Roll  Response 

-0.0109 

-0.01095 

for  comparing  the  results  of  the  numerical  linearization  with  Roskam's  tabulated  data 
is  to  form  plant  and  control  matrices  for  the  test  aircraft,  using  the  linear  algebraic 
method  taught  in  AE  3340  [Ref.  Sch  92].  The  resulting  A  and  B  matrices  were  then 
compared  by  applying  step  elevator,  rudder,  and  aileron  inputs  and  plotting  the  re- 
sults. The  results  for  the  step  elevator  input  are  given  in  Figure  4.13.  These  plots 
show  very  little  difference  in  the  vertical  velocity,  it;.  Differing  amplitudes  are  shown 
for  the  horizontal  velocity,  but  since  the  non-dimensional  velocity,  u/Vj,  from  the 
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dimensional  derivatives  was  scaled  to  be  equivalent  to  the  state,  u,  computed  in  the 
numerical  linearization,  these  magnitude  errors  are  not  indicative  of  a  poor  model, 
only  a  slight  difference  in  the  computed  damping  is  shown.  The  results  from  the 
analytic  model  were  scaled  up  by  the  nominal  airspeed,  Vj  to  compare  with  the  nu- 
merically linearized  results.  This  would  have  the  effect  of  magnifying  any  differences 
between  the  analytic  model  and  the  numerical  model.  The  natural  frequency  of  both 
the  analytic  and  numerical  results  are  quite  similar. 


40  60 

Time,  seconds 


Figure  4.13:  Comparison  of  Longitudinal  Responses  to  Step  Elevator  Input 

Lateral  responses  to  step  rudder  and  step  aileron  inputs  are  shown  in  Fig- 
ures 4.14  and  4.15.  Very  little  difference  between  the  models  is  visible  in  these  plots. 
The  errors  are  shown  in  Figures  4.16  -  4.1& 

It  can  be  concluded  that  the  results  from  the  numerical  linearization  are  quite 
close  and  furthermore  that  the  equations  used  in  the  model  are  indeed  correct.  More- 
over, the  linearized  results  presented  for  both  the  AROD  and  Bluebird  aircraft  should 
be  considered  accurate  and  are  suitable  for  the  Guidance,  Control,  and  Navigation 
system  designs  that  will  follow. 


52 


Plot  of  stop  Rudd©r  Data 

1 

200 

Ik 

V  nunnencal  

V  analytic 

p  numerical  — 
p  analytic 

rnumerxal 

r  analytic 

10O 
0 

\ 

- 

■100 

^^\^ 

- 

■200 

40  60 

Time,  seconds 


Figure  4.14:  Comparison  of  Lateral  Responses  to  Step  Rudder  Input 
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Figure  4.15:  Comparison  of  Lateral  Responses  to  Step  Aileron  Input 
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Figure  4.16:  Difference  in  Analytic  and  Numerical  Results,  Step  Elevator 
Input 


40  60 

Tune  seconds 


Figure  4.17:   Difference  in  Analytic  and  Numerical  Results,  Step  Rudder 
Input 
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Figure  4.18:  Difference  in  Analytic  and  Numerical  Results,  Step  Aileron 
Input 
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V.  SENSOR  AND  ACTUATOR  MODELING 

The  A  ROD  engine  and  actuators  were  modeled  as  a  second  order  transfer  func- 
tions, based  on  data  collected  by  Sandia  Labs  [Ref.  Wh  87].  Sensors  were  also  mod- 
eled as  second  order  transfer  function  based  on  data  supplied  by  Watson  Industries 
[Ref.  WAT  93], 

The  complete  non-linear  model  for  an  aircraft  should  include  models  of  the 
sensors  on  board  for  measurement  of  acceleration,  angular  rates,  pitch  and  roll  an- 
gles, and  headings.  The  inertial  device  chosen  for  the  AROD  project  was  the  Watson 
Industries  IMU-600D  inertial  measuring  unit.  This  device  contains  a  triaxial  ac- 
celerometer,  a  triaxial  rate  sensor,  two  liquid  pendulous  devices  (for  bank  and  pitch 
angle),  and  a  magnetic  heading  indicator.  The  characteristics  of  these  devices  must 
be  accurately  modeled,  since  it  is  the  sensor  output  that  drives  the  control  system. 
In  the  rest  of  this  section,  the  accelerometers,  rate  gyros,  and  inclinometers  will  be 
modeled,  as  well  as  the  sources  of  error  inherent  in  the  design  of  the  sensor  devices. 

A.       ACCELEROMETER  MODELING 

The  term  accelerometer  is  not  entirely  accurate,  since  the  device  does  not  mea- 
sure true  acceleration,  but  rather  the  difference  between  acceleration  and  gravity 
[Ref.  Bro  64].  This  effect  is  referred  to  as  the  Einstein  Uncertainty  Principle,  and  is 
represented  in  equation  form  as 

f  =  g-a  (5.1) 

where  g  and  a  are  the  specific  forces  of  gravity  and  acceleration  of  the  aircraft  and 
are  measured  positive  downwards.  The  accelerometer  model  relevant  to  this  equation 
is  shown  in  Figure  5.1.  The  tri-axial  accelerometer  of  the  IMU-600D  can  be  modeled 
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Figure  5.1:  Typical  Accelerometer  Model 

a^  three  simple  single-axis  accelerometers,  as  has  been  established  through  conversa- 
tions with  the  manufacturer  [Ref.  WAT  93].  A  schematic  representation  basic  device 
pictured  in  Figure  5.1  is  modeled  by  an  ordinary  differential  equation  [Ref.  Sil  91] 


..       /3  .       k 

X  -i X  H X  =  —y 

m         m 


(5.2) 


where  x  denotes  the  displacement  of  the  mass  from  its  equilibrium  position  and 
y  =  g  —  a  '\s  the  projection  along  the  case  axis  of  the  vector  sum  of  gravity  and 
acceleration.  The  terms,  /?,  A-,  and  m  represent  the  damping,  spring  coefficient,  and 
mass,  respectively,  of  the  device.  The  accelerometer  described  in  Equation  5.2  can 
be  modeled  as  a  second  order  low  pass  filter,  but  the  actual  accelerometer  has  a  flat 
response  up  to  lOOOHz,  so  it  was  not  modeled.  A  third  order  Chebychev  anti-aliasing 
filter  with  a  cut-off  frequency  of  20  Hz  was  added  to  the  accelerometers.  This  filter 
is  the  device  that  was  modeled.  The  Chebychev  filter  gives  the  advantage  of  a  flat 
passband,  and  a  very  sharp  drop  off  at  the  cut-off  frequency.  The  equation  used  to 
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describe  a  Chebychev  filter  [Ref.  St  88]  is 

\Hip{ju:)\'  = 


1 


(5.3) 


where  Cf^  is  the  Nth  order  Chebychev  polynomial,  e  is  the  parameter  thats  sets  the 
ripple  in  the  passband,  and  \Hip{juj)\^  is  the  magnitude  of  the  filter.  It  was  not 
necessary  to  compute  the  filter  analytically,  as  SiMULINK  provides  a  block  function 
which  performs  the  required  steps  based  on  the  passband  ripple  of  O.ldb  and  the 
desired  cutoff  frequency  of  20  Hz.  The  block  diagram  for  the  accelerometer  model 
is  shown  in  Figure  5.2.  Included  in  this  diagram  are  the  blocks  containing  the  error 
modeling,  as  well  as  the  blocks  used  to  correct  for  the  Einstein  uncertainty.  Figure  5.3 
shows  the  linear  synthesis  model  used  to  generate  the  accelerations  to  drive  the  sensor 
models.  The  synthesis  model  was  derived  from  the  Bluebird  model  discussed  in 
Chapter  IV. 
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Figure  5.2:  Accelerometer  Modeling 
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Figure  5.3:  Synthesis  Model  for  Accelerometers 

1.       Error  Model 

No  matter  how  well  the  sensor  device  is  constructed,  any  accelerometer 
is  subject  to  certain  errors  in  the  linear  acceleration  measurements.  These  errors 
can  occur  for  several  reasons;  some  as  mechanical  and  others  are  due  to  the  physical 
placement  of  the  accelerometer  on  the  aircraft.  The  mechanical  errors  accounted  for 
in  the  IMU-600D  tri-axial  accelerometers  are 

•  Acceleration  Bias 

-  average  readings  for  -\-lg  and  -l^^  loads 

•  Acceleration  Scale  Factor  error 

-average  difference  between  readings  from  +1^  and  -Ig  loads 


59 


•  Cross  Axis  Sensitivity  errors 

-measurements  due  to  the  misalignment  of  an  accelerometer  with  the  appropri- 
ate axis 
-measurements  of  off-axis  accelerations  are  measured 

•  Acceleration  Noise  Floor 

-threshold  below  which  measurements  of  acceleration  can  not  be  made 

Errors  in  the  measured  accelerations  can  also  occur  if  the  accelerometers 
are  not  located  at  the  e.g.  of  the  aircraft,  since  arbitrarily  located  points  on  a  body 
will  experience  additional  accelerations  due  the  the  angular  momentum  of  the  body. 
A  mathematical  correction  for  the  error  will  be  presented  later,  and  will  describe  the 
angular  and  linear  accelerations  of  an  arbitrary  point  on  a  rigid  body.  However,  the 
correction  will  not  be  applied  to  the  models  here,  since  the  correction  is  expected  to 
have  a  very  small  effect  on  the  sensor  measurements  due  to  the  small  displacements 
away  from  the  aircraft  e.g. 

Error  terms  are  quantified  in  terms  of  full-scale  measurements.  The  me- 
chanical errors  are  tabulated  in  Table  5.1  along  with  other  important  characteristics. 
Figure  5.4  shows  the  block  diagram  with  error  inputs  applied  to  the  accelerations 

TABLE  5.1:  ACCELEROMETER  CHARACTERISTICS 


Acceleration  Range 

±2g's 

Acceleration  Bandwidth 

20  Hz 

Acceleration  Bias 

0.2%  of  Full  Scale 

Acceleration  Scale  Factor 

0.2%  of  Full  Scale 

Acceleration  Noise  Floor 

0.0005  g's 

Cross  Axis  Sensitivity 

0.5%  of  Full  Scale 
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Figure  5.4:  Error  Model  for  Accelerometers 

measured  by  the  accelerometers.  The  output  from  the  error  computations  in  Fig- 
ure 5.4  is  the  measured  acceleration  from  the  accelerometer  output  to  the  control 
system.  The  cross  axis  terms  are  determined  through  a  matrix 


X  = 


0     Cy    e, 

Cx       0       Cj 

Cx         Cy  0 


(5.4) 


where  e,  is  the  cross  axis  error  term  and  x  is  the  error  in  the  acceleration  due  to  the 
cross  axis  sensitivity. 

2.       Results  and  Validation 

Accelerations  were  measured  for  step  aileron,  elevator,  and  rudder  inputs, 
and  the  results  then  compared  to  the  actual  accelerations  computed  for  the  equations 
of  motion  for  the  aircraft.  A  linear  synthesis  model  was  used  for  the  initial  testing, 
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while  the  results  from  the  non-linear  model  were  used  for  validation  of  the  accelerom- 
eter  model.  Figures  5.5,  5.6,  and  5.7  show  comparisons  between  measured  and  actual 
accelerations  generated  from  simulations  of  the  non-linear  model.  The  accelerations 
computed  for  the  longitudinal  and  lateral  cases  were  in  close  agreement  with  the 
accelerations  computed  by  the  non-linear  model. 
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Figure  5.5:  Measured  and  True  Acceleration  From  a  Step  Elevator  Input 
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Figure  5.6:  Measured  and  True  Acceleration  From  a  Step  Aileron  Input 
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Figure  5.7:  Measured  and  True  Acceleration  From  a  Step  Rudder  Input 
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B.     RATE  GYRO  MODELING 

The  rate  gyros  consist  of  a  rotating  disk,  mounted  in  a  gimbal  mechanism. 
Though  both  single  and  two  degree  of  freedom  gyros  are  common,  the  tri-axial  angle 
rate  gyro  supplied  in  the  IMU-600D  is  modeled  as  three  single-degree-of-freedom 
rate  gyros,  each  measuring  the  angular  rate  along  a  particular  axis.  The  dynamics  of 
a  gyro  can  be  modeled  [Ref.  Sil  91]  using  Euler's  law 

^Ng  =^  Lg 

wher  {B]  and  {G}  are  the  body  and  gyro  coordinate  systems,  respectively.  This  can 
be  expanded  to 

^Ng  =  ^u;g  X  ^Lg  +  g^l^^G 

where  '^Lg  =  Ig^^g  and  the  time  derivative,  4:^ Lg  is  zero  when  the  wheel  rotates 
at  a  constant  speed.  Thus,  except  for  the  period  when  the  wheel  is  coming  up  to 
speed,  the  equation  for  gyroscopic  motion  in  a  rate  gyro  can  be  written  as 

^Ng  =  ^u:g  X  ^Lg  (5.5) 

where  ^Ng  is  the  torque  vector  acting  on  the  gyro  element  and  ^uJg  is  the  angular 
rate  of  the  gyro  frame.  Transfer  functions  can  be  developed  for  the  torque  input  to 
the  input  axis  as  shown  in  Figure  5.8  and  the  rate  output  of  the  output  axis.  This 
derivation  was  not  developed  since  data  required  for  the  gyro  disk  inertia,  Ig  and  the 
speed  of  rotation,  ^wg,  among  other  terms,  is  not  available  from  the  manufacturer. 
The  rate  gyros  were  modeled  as  second  order  transfer  functions,  with  u,'n  =  50  Hz. 
A  third  order  Chebychev  filter  with  a  cut-off  frequency  of  20  Hz  was  added  to  the 
gyro  model  as  an  anti-aliasing  filter.  This  Chebyshev  filter  was  identical  to  the  ones 
described  for  the  accelerometers.  The  block  diagram  of  the  rate  gyro  model  is  shown 
in  Figure  5.9.    The  linear  synthesis  model  is  shown  in  Figure  5.10  and  was  used  to 

64 


Magnetic 
suspension' 


Signal 
generator 


Gyroscope 
element 


Torque 
glnwator 


Magnetic 
suspension 

Case 
(gyro  housing) 


Output  axis 
(O) 


Spin  9»s(s) 


angular  ^■ 
momentunr» 
<H) 


Figure  5.8:  Functional  Diagram  of  a  Rate  Gyro  [Ref.  Bro  64] 

test  the  gyro  models.    In  Figure  5.9  the  gyro  elements  are  shown  for  each  axis  and 
blocks  for  the  error  calculations  are  shown  as  well. 
1.       Error  Modeling 

Error  terms  are  also  present  in  the  rate  gyros  and  are  due  to  either  physical 
location  on  the  aircraft,  or  mechanical  errors,  as  in  the  accelerometers.  Errors  due 
to  physical  placement  away  from  the  e.g.  can  be  corrected  by  using  the  equations 
derived  for  the  linear  and  angular  acceleration  of  an  arbitrary  point  on  the  aircraft, 
shown  later  in  this  chapter.  Mechanical  errors  are  defined  in  the  same  manner  as  was 
done  for  accelerometers  in  Section  A,  and  are  listed  in  Table  5.2  along  with  other 
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Figure  5.9:  Rate  Gyro  Model 


TABLE  5.2:  GYRO  CHARACTERISTICS 


Rotational  Rate  Range 

±\UMtglstc 

Rotational  Rate  Bandwidth 

20  Hz 

Rotational  Rate  Bias 

2.0%  of  Full  Scale 

Rotational  Rate  Scale  Factor 

0.5%  of  Full  Scale 

Rotational  Rate  Noise  Floor 

0.05%  of  Full  Scale 

Cross  Axis  Sensitivity 

0.5%  of  Full  Scale 

important  characteristics.  The  cross  axis  error  is  modeled  in  the  same  way  as  for 
the  accelerometers  by  the  use  of  the  same  matrix  for  computing  errors  in  the  angular 
rates. 

2.       Results  and  Validation 

Angular  rates  were  measured  for  step  aileron,  elevator,  and  rudder  inputs 
and  compared  to  the  actual  angular  rates  computed  from  the  equations  of  motion. 
First  a  linearized  synthesis  model  was  used  in  the  testing  stages;  then  the  sensors 
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Figure  5.10:  Rate  Gyro  Synthesis  Model 

were  integrated  with  the  non-linear  aircraft  simulation.  Comparisons  of  actual  and 
measured  acceleration  are  given  in  Figures  5.11,  5.12,  and  5.13.  These  figures  show 
the  accuracy  of  the  angular  rate  sensors  in  measuring  the  angular  rates  computed 
with  the  non-linear  equation  of  motion  model. 
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Figure  5.11:  Measured  and  True  Angular  Rates  From  a  Step  Aileron  Input 
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Figure  5.12:    Measured  and  True  Angular  Rates  From  a  Step  Elevator 
Input 
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Figure  5.13:  Measured  and  True  Angular  Rates  From  a  Step  Rudder  Input 
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C.     PITCH,  ROLL,  AND  HEADING  SENSOR  MODELING 

The  last  group  of  sensors  to  be  modeled  are  those  which  measure  the  pitch,  roll, 
and  heading  angles.  The  pitch  and  roll  angle  measurements  are  made  with  liquid 
pendulous  devices.  These  are  devices  that  have  an  electrolyte  contained  in  a  vial, 
and  by  measuring  the  capacitance  changes  of  the  vial  as  the  electrolyte  moves  in 
response  to  aircraft  angle  changes,  the  pitch  and  roll  angles  can  be  determined.  The 
heading  sensor  is  a  magnetic  heading  sensor. 

The  primary  concern  with  angular  measurements  is  that  accurate  readings  are 
obtained,  even  when  the  aircraft  is  experiencing  linear  and  angular  accelerations. 
When  these  errors  cannot  be  corrected,  the  control  system  must  be  able  to  compen- 
sate for  the  errors.  An  inclinometer  is  typically  modeled  as  a  pendulum,  shown  in 
Figure  5.14  attached  to  a  block,  which  can  be  considered  to  represent  an  aircraft. 
The  equations  describing  the  motion  of  the  pendulum  are  derived  in  detail  later  in 
this  section.  The  transfer  function  for  the  pendulum  inclinometer  can  be  represented 
as  a  second  order  transfer  function  with  ^'n  —  0.8  Hz  and  (,'  =  0.5.  [Ref.  W.^T  93] 

0  5.03^ 

0  ~  52 +  50.35 +  5.032" 

1.       Error  Modeling 

The  inclinometers  are  subject  to  several  sources  of  errors,  from  both  linear 
and  angular  accelerations,  and  from  mechanical  imperfections.  Mechanical  errors  are 
listed  in  Table  5.3,  and  modeled  as  shown  in  Figure  5.15.  Errors  due  to  angular 
velocity  can  be  compensated  for  with  a  complementary  filter.  Selecting  the  time  con- 
stants appropriately  will  allow  direct  measurements  of  the  angle  from  the  inclinometer 
at  low  frequencies,  and  from  integrated  angular  rates  at  higher  frequencies  where  the 
inclinometer  is  inaccurate. 
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Figure  5.14:  Simple  Pendulum  Inclinometer 

TABLE  5.3:    INCLINOMETER  AND  HEADING  SENSOR  CHARAC- 
TERISTICS 


Pitch  and  Roll  Range 

±50  deg 

Pitch  and  Roll  Bandwidth 

1/2  Hz 

Pitch  and  Roll  Accuracy 

0.2  deg 

Heading  Range 

±180  deg 

Heading  Accuracy 

3.0  deg 

Heading  Repeatability 

0.5  deg 

Heading  Linearity 

0.5% 

It  was  determined  that  the  effects  of  linear  acceleration  were  also  important 
to  model.  In  order  to  model  these  effects  it  was  necessary  to  derive  a  transfer  function 
from  the  aircraft's  linear  acceleration  to  the  angle  caused  by  that  acceleration.  The 
starting  point  was  to  model  the  inclinometer  as  shown  in  Figure  5.14.  The  derivation 
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Figure  5.15:  Inclinometer  Error  Modeling 

results  in  two  coupled  equations,  the  first  of  which  will  be  ignored,  since  the  equations 
of  motion  for  the  aircraft  are  known  and  the  motion  of  the  inclinometer  should  have 
little  effect  on  the  aircraft  motion.  The  second  equation  is  for  the  motion  of  the 
pendulum  as  influenced  by  the  aircraft's  acceleration  and  is  what  was  used  for  the 
linear  acceleration  errors. 

First  it  was  necessary  to  define  the  coordinates  used  to  describe  the  incli- 
nometer system  as  g  =  {x,6)  (see  Figure  5.14).  Next,  since  Lagrangian  methods  were 
used  for  the  derivation,  the  kinetic  energy,  T,  of  the  system  was  defined  as 


T    =     l/2Mi^  + l/2m(i-  +  /^cos^)^- l/2m(/^sin^)^ 
=    1/2MP  +  1/2j7}P^  +  "^  I '2m  Pi, 


(5.6) 


where  Pj.  and  Py  represent  the  position  of  the  pendulum  in  Cartesian  coordinates. 
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The  potential  energy,  V,  of  the  pendulum  can  be  written  as 

V  =  V^  =  mgl{\  -  cos6>),  (5.7) 

where  M,  m,  and  1  are  the  mass  of  the  aircraft,  mass  of  the  pendulum,  and  distance 
from  the  aircraft  e.g.  to  the  pendulous  mass.  The  terms, P,  and  6,  represent  the 
position  vector  to  the  pendulous  mass  and  the  angle  made  by  the  pendulum  with  the 
vertical  plane.  The  governing  equation  for  Lagrangian  dynamics  is 

i^_^^g,  (5.8) 

dt  dq,      dqi 

where  C  is  the  Lagrangian  operator,  C  =  T  —  V^  and  Q,  represents  the  non- 
conservative  forces  acting  on  the  body.  After  some  rearranging, 

C  =  1/2(M  +  m)i'^  +  mlid cosd  +  l/2ml^e^  -  mgl{\  -  cosO)  (5.9) 

The  result  in  Equation  5.9  can  be  substituted  into  Equation  5.8  resulting  in 

F)  r 

-—    =    {M -^  m)x  +  ml6 cos6 
ox 

— -     =    mlxcosO  -\- 1^9 
06 

—    -    0 
dx 

dC 

-r-r     =     —mlx9  s\n6  —  mqlsmO 
do 

rl   Pi  r 

—  -—    =    (M -^m)x -{-mlcosOO  -  mls\n9  9^ 
dt  ox 

d  dC  ..  .       •  2" 

-; =    mlcosOx  —  mlsinOO -\- ml  9. 

dtdO 

When  these  partial  derivatives  are  substituted  into  Equation  5.8,  the  resulting  equa- 
tions of  motion  for  the  pendulum  are 

(M  +  m);r  +  m/cos^^  -  m/sin^^^     =     Qi 
ml cos9x  +  mP9 -\- mglsin9    =    Q21 
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where  Qi  and  Q2  are  terms  representing  damping  in  the  system.  The  term  Qi  repre- 
sents damping  on  the  aircraft  for  Qi  =  —fi^i  and  the  term  Q2  represents  the  viscous 
damping  of  the  pendulous  mass,  Q2  =  —f3r0.  The  equations  can  be  linearized,  where 

{M -\- m)x -\- mW -\- 0rx    =    0  (5.10) 

mlx-\-mPe-\-  l3re-\-mgW    =    0.  (5.11) 

Now  it  is  apparent  that  Equation  5.10  is  completely  determined  by  the  aircraft  equa- 
tions of  motion  that  have  already  been  modeled.  A  Laplace  transform  of  Equa- 
tion 5.11  can  be  performed  to  find  the  desired  transfer  function 

g=  -'"  .  (5.12) 

where  the  similarity  to  the  standard  second  order  transfer  function  is  noted,  making 
the  term  g/l  =  lo^.  The  term  —1//  in  the  numerator  can  be  solved  for  by  substituting 
g  =  32.174  f/s^  and  LOn  =  5.027  rad/sec,  or  O.SHz,  with  the  result  /  =  1.27  /.  This 
result  is  then  substituted  into  the  block  diagram  for  the  inclinometer  models  shown 
in  Figure  5.15.  Responses  for  step  inputs  are  shown  in  Figures  5.16,  5.17,  and  5.18. 
In  the  lateral  cases,  there  is  very  little  difference  between  the  measured  and  actual 
angles.  The  longitudinal  case  is  quite  different.  In  Figure  5.16.  there  is  a  considerable 
difference  between  the  actual  pitch  angle,  0,  and  the  pitch  angle  measured  with  the 
inclinometer  model.  This  difference  the  linear  acceleration  of  the  aircraft  as  a  step 
elevator  input  is  applied,  and  must  be  compensated  for  before  a  reliable  control  system 
can  be  developed. 
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Figure  5.16:  Response  to  a  Step  Elevator  Input 
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Figure  5.17:  Response  to  a  Step  Rudder  Input 
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Figure  5.18:  Response  to  Step  Aileron  Input 


The  error  in  inclinometer  readings  due  to  linear  acceleration  has  been  de- 
termined experimentally  by  the  manufacturer,  and  can  be  represented  as  a  second 
order  transfer  function  from  g's  to  0, 

e  ^  157.08^ 

g        ^2.^  157.08s +  157.082'  ^  "     ^ 

where  the  gain,  A',  was  determined  from  an  input  of  1.2  V  at  10  mV/^,  and  the 
resulting  output  of  0.4  V  at  60  mV/deg.  Thus  the  gain.  A',  was  17.99. 

D.     MODELING  OF  AN  ARBITRARY  SENSOR  PLACEMENT 

The  actual  sensor  placement  in  either  the  Bluebird  or  the  AROD  aircraft  is  not 
at  the  e.g.,  as  was  assumed  in  the  previous  sections.  In  order  to  model  the  actual  linear 
and  angular  accelerations  at  the  sensors,  regardless  of  the  position  on  the  aircraft, 
new  equations  of  motion  must  be  derived.  These  equations  can  then  be  applied  to 
the  sensor  inputs  to  obtain  a  proper  output  from  the  sensor.  First  the  equations 
for  linear  acceleration  will  be  derived  using  Newton's  third  law  for  conservation  of 
momentum,  then  the  equations  for  angular  acceleration  will  be  derived  using  Euler's 
law  for  conservation  of  angular  momentum.  These  equations  must  be  expressed  in 
terms  oi  [B]  since  all  the  information  measured  by  the  sensors  will  be  determined  in 
the  body  coordinate  system. 

1.      Linear  Accelerations 

In  the  derivation  for  equations  of  motion  of  an  arbitrary  point  on  the  air- 
craft, the  coordinate  systems  representing  the  inertial  and  body  coordinate  systems, 
similar  to  those  shown  in  Figure  3.1  are  used.  Supplementing  the  derivation  of  equa- 
tions of  motion  for  an  aircraft,  the  motion  of  the  sensor  location  on  the  aircraft,  given 
by  Pq  is  necessary.  In  the  inertial  coordinate  system,  [U],  the  position  of  Q  can 
be  written  as  ^^ Pq  —  ^R^Pq  -\-  ^ Pbo-   The  velocity  is  first  determined  by  applying 
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Coriolis'  theorem  from  Equation  3.24  to  obtain  the  first  derivative.  Here, 

" Vb  S  "P,  =  "Pbo  +  "jfiR^PQ),  (5.14) 


where 


U  "  (U  dBd    \  —  B  "■  iU  dBd    \    I    V    .      ,,  ,U  dB 


j^CbR"Pq)  =     J^R^Pq)  +  '^B  X  (^i?^Pg).  (5.15) 


The  terms  ^'^  and  ^^  refer  to  derivatives  taken  with  respect  to  the  inertial  and  body 
coordinate  systems,  respectively.  The  resulting  expression  for  velocity  is 

%  =  ^Vbo  +  ^x(^5/?%),  (5.16) 

where  ^^(gi^^Pg)  =  0  since  Q  is  a  point  fixed  on  the  aircraft.    Accelerations  are 
derived  by  applying  Equation  3.24  twice  to  Equation  5.16 

''Vq  =  ^j/'^Bo  +  %  X  {%R%))  +  %  X  {"Vbo  +  %  X  (^i?%)).  (5.17) 

The  ^^(O  term  is  differentiated,  resulting  in 

^^(^'Vso  +  %  X  {U^Pq})  =  ^J/Vbo)  +  "d;B  X  ^it>%  +  ^  x  ^^Fg.        (5.18) 

Substituting  the  results  from  Equation  5.18  into  Equation  5.17,  the  result  is 

""Vq  =  ^j/^Bo)  +  2  ''u:b  X  "Vq  +  ''iCB  X  ^/?%  +  ^  x  (^  x  %R^Pq).       (5.19) 

The  desired  result  in  {B}  is  obtained  by  simply  premultiplying  Equation  5.19  by  gR 

^Vq  =  ^l^VBo)  +  2  ^u;b  X  ^Vq  +  ^UB  x  %  +  «c.'b  x  (^..-b  x  %),         (5.20) 
at 

where  the  identity  [Ref.  Sp  89] 

^R{lob  X  ^Vq)  =  (^ilc^'fi)  x  {f,R"VQ) 

was  used.    This  result  can  now  be  substituted  into  Equation  3.27  and  the  resulting 
expression  can  be  solved  for  -jii^VBo)- 
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2.      Angular  Accelerations 

No  further  work  is  needed  to  derive  the  expressions  for  angular  velocity 
and  acceleration  at  a  arbitrary  point,  since  for  a  rigid  body,  these  quantities  remain 
the  same  anywhere  on  the  body. 

iOBO  =  u;q^Q  e  {B] 
since  a;  x  u;  =  0. 


79 


VI.  CONCLUSIONS  AND 
RECOMMENDATIONS 

A.  CONCLUSIONS 

Based  on  the  data  presented  in  this  thesis,  the  following  conclusions  are  made: 

•  High-fidelity  models  of  both  the  AROD  and  Bluebird  aircraft  were  implemented 
using  SiMULINK  software.  Use  of  the  block  diagram  structure  of  the  model 
allows  changes  to  be  easily  made  and  requires  no  programming  ability,  other 
than  the  use  of  Matlab. 

•  These  models  accurately  represent  the  sensors,  actuators,  and  engines  associated 
with  the  particular  aircraft. 

•  Expected  errors  in  the  accelerometers  and  rate  gyros  are  very  small.  Pitch 
errors  from  the  inclinometers,  due  to  linear  accelerations,  will  be  much  greater. 

•  The  models  that  were  developed  only  represent  one  flight  condition.  The  AROD 
was  modeled  only  in  a  hover  and  the  Bluebird  was  only  modeled  in  a  cruise 
condition.  Further  work  on  the  control,  guidance,  and  navigation  systems  will 
be  concentrated  in  these  flight  regimes.  The  data  files  for  a  given  flight  condition 
are  easily  replaced  with  tables,  when  more  test  data  are  available. 

B.  RECOMMENDATIONS 

Based  on  the  conclusions  presented  above,  and  the  problems  encountered  during 
this  study,  the  following  recommendations  are  made. 

•  The  Department  of  Aeronautics  and  Astronautics  should  update  the  UNIX  labs 
to  include  SiMULINK  with  Matlab  4.0.  This  will  allow  increased  instructional 
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use  in  flight  dynamics,  controls,  and  avionics  courses.  Additionally,  a  practical 
research  tool  would  be  available. 

•  The  work  in  this  thesis  must  be  modified  to  include  the  aerodynamic  charac- 
teristics of  the  Archytas  aircraft.  Achieving  this  step  will  require  wind  tunnel 
testing  for  the  Archytas  in  all  the  expected  phases  of  flight,  especially  the  ex- 
tremely non-linear  transition  phase  from  vertical  to  horizontal  flight. 

•  Integration  of  the  quaternion  based  rotation  matrix  should  be  completed  in 
order  to  take  advantage  of  the  superior  qualities  of  quaternions  in  computational 
models. 
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APPENDIX  A:  MATHEMATICAL 
PROPERTIES 

In  this  appendix,  some  of  the  mathematical  properties  used  in  the  text  are 
described. 

A.     CROSS  PRODUCT  PROPERTIES 

In  this  section,  important  properties  of  the  cross  product  and  cross  product 
matrices  are  described  The  cross  product  between  two  vectors, 


A  = 


'   Gr   ' 

«!/ 

.    ^^    . 

and  B 


is  defined  by 


Ax  B  = 


dyb^  —  a^by 

Qj-by      —     ttybj. 


Properties  of  the  cross  product  are: 


A  X  B  =  {Ax)B,  where  Ax  = 
and  is  called  a  cross  product  matrix 


■      0 

—  G; 

Qy 

a. 

0 

-Ox 

.   -ay 

flx 

0 

•  iR{V  X  U)  =  iiRV)  X  {-^RU)  if  iR  is  a  rotation  matrix  with  V  and  U  in  the 
same  coordinate  system.  That  this  matrix  distributes  across  the  cross  product 
is  obvious  since  rotation  matrices  preserve  space  geometry. 

•  B^i^'  '  ^)  —  [b^^')   '   {'b^^')  ^^^  ^he  same  reasons  as  above. 


Ax{B  xC)  =  {A-C)B-{A-B)C 
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•  A  X  {B  X  C)  =  B  X  {A  X  C)  +  C  X  {B  X  A) 

•  Ax  B  =  -B  X  A 

•  —Ax  =  A^x 

B.     DERIVATIVES  OF  VECTORS 

For  any  free  vector,  ^Vq  (i.e.  a  velocity  vector  and  any  rotation  matrix,  ^i?,  the 
derivative  of  the  velocity  of  Q  computed  in  {B}  and  expressed  in  {^1}  and  denoted 
as  ^{^Vq)  is  given  as  [Ref.  Sil  91] 

since  as  shown  in  [Ref.  Cra  86], 

SO  then 

iR^^R  =  ^nBX 
The  same  process  can  be  carried  out  for  an  angular  velocity  vector  "^Qb-^ 

=    iRf/i'QB)) 

And  if  the  origins  of  the  coordinate  systems  {A]  and  {B]  are  coincident,  the  derivative 
can  be  expressed  as 

~  'd'v 
83 


For  an  applied  vector  ^Pq  (i.e.  the  position  of  point  Q  in  the  {B}  coordinate  system) 
and  a  rotation  matrix  gR,  the  time  derivative  of  the  position  vector  of  Q,  expressed 
in  {A}  ,  ^("^-Pq)  as  a  function  of  it's  derivative  computed  in  {B]  is  given  by  [Ref. 
Sil  91] 

=    ^BR%  +  BRf^{^PQ)  +  ^yB 

Therefore  the  velocity  of  the  point  Q  can  be  expressed  in  {A}  as 

^Vq  =  ^^B  X  {JR  %)  +  iR  ^Vq  +  Hb 
or  expressing  the  velocity  of  Q  in  the  {B]  coordinate  system 

In  the  case  where  the  origins  of  {A]  and  {B]  are  coincident  then  the  resulting  ex- 
pression is 

^VQ  =  ^^B  X  %  +  ^Vq  +  ^VB. 

C.     EQUATIONS  USED  FOR  LINEARIZATION 

For  any  vector  equation,  H{c)  =  a  x  b  where  c  =  [a  b]  and  a,  b,  and  c  are  all 
vectors,  the  Taylor  series  approximation  can  be  written  as 

dH 

H{c)  =  H{co)  ^  ^\c=cSc^  H.O.T. 
oc 


and 


dH       dH^        ,        dH.       ,, 


Now  dH/da  can  be  written  as 


dH       d{a  X  6)       d(-b  x  a)       di-b)  ,       da 

-^-  =  — ^ =  ^ =  -T. —  X  a  -  6  X  —  =  -6x 

oa  oa  oa  oa  da 
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Using  the  same  relations,  dH/db  can  be  written 

dH       d{-bxa) 


db  db         ^^" 
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APPENDIX  B:    NUMERICAL  RESULTS 


A.     AROD  RESULTS 

The  following  results  are  from  the  numerical  linearization  of  the  kinematic  equa- 
tions of  motion  for  the  AROD.  Using  the  trim  command,  based  on  an  initial  ©o  =  7r/2, 
resulted  in  the  state  vector  and  input  vector  of: 


X  = 


0 
0 
0 
0 
0 
0 
0 
1.5700 
0 


and  u  = 


The  linmod  command,  [a,b,c,d]:=linmod('baslnmod\x,u),  produced  the  following  lin- 
earized system: 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-1.5024 

0 

0 

0 

0 

1.5112 

0 

b  = 


10  0  0  0  0 
0  10  0  0  0 
0  0  10  0  0 
0  0  0  10  0 
0  0  0  0  10 
0    0    0    0    0    1 


and  the  c  and  d  matrices  were  empty,  since  no  outputs  were  defined. 

The  following  results  are  for  the  numerical  linearization  with  gravitational  ef- 
fects added  to  the  2-3-1   Euler  angle  kinematic  equation  model.    The  trim  com- 
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mand,  using  the  same  initial  conditions  as  before,  resulted  in  the  following  vectors: 
[x,u]=trim('gravImod',xO,[],0,ix,[],0) 


X  = 


-0.00 
0.00 

-0.00 
0.00 
0.00 
0.00 
0.00 

1.5700 
0.00 


and  u  = 


32.1740 
-0.00 

-0.0256 

-0.00 

0.00 

-0.00 


Linearizing  the  system,  based  on  the  state  and  input  vector  found  by  trimming  at 
the  desired  conditions  resulted  in:  [a,b,c,d]=linmod('gravlmod\x,u) 


0 

0.00 

0     0   0.00 

0.00     0 

-0.0255 

0.0002 

-0.00 

0  0.00  -0.00     0 

0.00  0.0256 

-0.00 

32.1740 

0.00 

-0.00 

0  -0.00  -0.00 

0  -0.00 

-32.1740 

0 

0 

0 

0     0  -0.00 

-0.00     0 

0 

0 

a  = 

0 

0 

0   0.00     0 

-1.5204     0 

0 

0 

0 

0 

0  -0.00  1.5112 

0     0 

0 

0 

0 

0 

0   1.00  -0.00 

0.00     0 

0 

-0.00 

0 

0 

0     0   1.00 

-0.00  -0.00 

0 

0.00 

0 

0 

0     0   0.00 

1.00  -0.00 

0 

0 

and 

6  = 

■  1.00    0    0 
0  1.00    0 
0    0  1.00 

0      0      0  ■ 

0    0    0 
0    0    0 

0    0    0 

1.00    0    0 

' 

[ 

0    0    0 

0  1.00    0 

0    0    0 

0    0  1.00 

Trimming  the  full  EOM  model  at  hover  conditions  resulted  in  the  following: 

-0.00 
-0.00 

0.00 
-0.00 
X  =  0.00       and  u  = 

-0.00 

0.00 
1.5708 

0.00 


0.00 

0.00 

0.1807 

6387.2 
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The  subsequent  linearization  for  the  complete  model  yielded: 


a  — 


and 


0 

0 

0 

0 

0 

0 

0 

0.0002 

0.0002 

0.00 

0 

0 

0.00 

0 

0.00 

0.00 

0.00 

32.1740 

0 

0 

0 

0.00 

0 

0 

-0.00 

-32. 

1740 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-0.00 

0 

-1.5174 

0 

0 

0 

0 

0 

0 

-0.00 

1.5082 

0 

0 

0 

0 

0 

0 

0 

1.00 

— 

o.oo 

0 

-0.00 

0 

0 

0 

0 

0 

0 

1.00 

— 

0.00 

0.00 

0 

0.00 

0 

0 

0 

0 

0 

1.00 

0.00 

0 

0 

- 

0 

0 

0 

0.0112  ] 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

24.8194 

0.0003 

h 

zr 

-6.6192 

0 

0 

0.00 

0 

-6.5791 

0 

0.00 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0  . 

Computation  of  the  eigenvalues  gave: 


eig(a)  = 


0 
0 
0 
0 
0 
-0.00 

0+  1.5128?; 

0-1.5128? 
0 


Now  consider  the  quaternion  based  model  Rpm  is  fixed  at  6800  RPM. 

To  hold  the  desired  vertical  attitude,  the  quaternion  states  were  held  fixed: 

7 
8 
9 


IX  = 


10 
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where  the  state  vector  of  nominal  conditions  was: 


2-0  = 


0 
0 
0 
0 
0 
0 

0.7071 
0 

0.7071 
0 


Trim  and  linearization  of  the  quaternion-based  gravitational  model  resulted  in  the 
following  state  and  input  vectors,  as  well  as  linearized  a  and  b  matrices: 


X  = 


-0.00  ■ 

0.00 

-0.00 

■  32.1740 

0.00 

-0.00 

-0.00 
0.00 

,  u  = 

-0.00 
0.00 

0.7071 

0.00 

-0.00 

0.00 

0.7071 

0.00 

a  =,  columns  1-6, 


0 

0.00 

0.00 

0 

0.00 

0.00 

-0.00 

0 

0.00 

-0.00 

0 

0.00 

-0.00 

-0.00 

0 

-0.00 

-0.00 

0 

0 

0 

0 

0 

0.00 

-0.00 

0 

0 

0 

0.00 

0 

-1.6051 

0 

0 

0 

0.00 

1.6062 

0 

0 

0 

0 

0.00 

-0.3536 

-0.00 

0 

0 

0 

0.3536 

-0.00 

0.3536 

0 

0 

0 

0.00 

0.3536 

0.00 

0 

0 

0 

-0.3536 

-0.00 

0.3536 
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columns  7-10 


-45.5009 

0.00 

45.5012 

0 

0 

0 

0 

0.00 

-0.00 

0.00 


0 

-45.5009 

-0.0003 

0 

0 

0 

-0.00 

0 

-0.00 

-0.00 


-45.5009 

0.00 

-45.5012 

0 

0 

0 

0.00 

0.00 

0 

-0.00 


0 

45.5009 

0.0003 

0 

0 

0 

-0.00 

0.00 

0.00 

0 


and 


b  = 


1.00  0         0 

0  1.00         0 

0  0  1.00 

0  0         0    1.00 

0  0         0 

0  0         0 


0 

0 

0  " 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1.00 

0 

0 

0 

1.00  . 

Eigenvalues  are  computed  as: 


eig(a) 


-0.00  +  O.OOi 

-0.00  -  0.00^ 

0.00 

0 

-0.00  +  1.60572 

-0.00-  1.60572 

-0.00 

0.00 

-0.00  +  0.00? 

-0.00  -  0.00?' 


B.     BLUEBIRD  NUMERICAL  RESULTS 

Results  for  the  numerical  trim  and  linearization  of  the  Bluebird  model  are  given 
below.    The  results  are  from  the  trim  of  the  kinematic  equations  of  motion  for  the 
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Bluebird  model  are  based  on  the  nominal  conditions  of: 

72.9954 

0 

6.6757 

0 

xO=  0 

0 

0 

0.0912 

0 

with  the  states,  u,  ty,  and  0  held  constant,  ix  =   [15  8]'.    The  trim  command, 

[x,u]=trim('basic',xO,[],[],ix, [],[])  resulted  in  the  following  state  and  input  vectors: 


X  = 


72.9954 
0 

6.6757 
0 
0 
0 
0 

0.0912 
0 


and  u  = 


The  numerical  linearization  of  the  basic  model,  [A,B,C,D]=linmod('basic\x,u),  re- 
sulted in  the  following  A  and  B  matrices: 

0    0    0  0    -6.6757  0 

0    0    0    6.6757  0    -72.9954 

n    n    n  n      70  oqc 


0    0    0 

0 

72.9954 

0 

0    0    0 

0 

0 

0 

0    0    0 

0 

0 

0 

0    0    0 

0 

0 

0 

and 

10    0    0    0    0 

0  10    0    0    0 

0  0    10    0    0 

0  0    0    10    0 

0  0    0    0    10 

0  0    0    0    0    1 

The  C  and  D  matrices  were  empty  since  no  outputs  were  defined. 


B  = 
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Results  for  trim  and  linearization  of  the  kinematic  model  with  gravitational 
effects  added  yielded  the  following  state  and  input  vectors: 


X  = 


72.9954 

0.00 

6.6763 

0.00 

0.00 

-0.00 

-0.00 

0.0912 

0.00 


and  u  = 


2.8772 

0.00 

-31.4605 

-0.00 

0.3945 

-0.00 


The  numerical  linearization  of  the  model  resulted  in  the  following  A  and  B  matrices: 
A  = 


-0.0007 

-0.00 

0.0079 

0 

-6.6763 

0.00 

-0.00 

-32.0451 

0 

0.00 

0 

0.00 

6.6763 

0 

-72.9954 

32.0403 

0.00 

0 

-0.0001 

0 

0.0007 

-0.00 

72.9954 

0 

-0.0002 

-2.8773 

0 

-0.00 

0.0087 

-0.00 

0 

-0.00 

0.00 

-0.00 

-0.00 

0 

-0.00 

0.00 

0.0005 

0.00 

0 

-0.00 

0.00 

0.0361 

0 

-0.00 

0.0010 

0.00 

-0.00 

-0.00 

0 

-0.00 

-0.00 

0 

0 

0 

0 

1.00 

-0.00 

0.0915 

0.00 

-0.00 

0 

0 

0 

0 

0 

1.00 

0.00 

0.00 

0 

0 

0 

0 

0 

0 

-0.00 

1.0042 

0.00 

-0.00 

0 

and 

10    0    0    0  0 

0    10    0    0  0 

0    0    10    0  0 

0    0    0    10  0 

0    0    0    0    10 

0    0    0    0    0  1 

Results  for  trimming  the  entire  model  were  A  = 


B 


■0.0622 

0.00 

-0.7558 

0.00 

0.0155 

0.00 

0 

0 

0 


0.00 

-0.3865 

0.00 

-0.1455 

0.00 

0.1418 

0 

0 

0 


0.3431 

0.00 

-4.7145 

0.00 

-0.1907 

0.00 

0 

0 

0 


0 

1.7450 

0.00 

■5.3700 

0.00 

-1.0589 

1.00 

0 

0 


-1.6187 
0 

67.1233 
0.00 

-3.1266 
0.00 
0.00 
1.00 
0.00 


0.00 

-71.6817 

0 

1.5003 

0.00 

-0.7970 

0.0915 

0.00 

1.0042 


0.00 
32.0403 
-0.0002 
0.00 
0.00 
0.00 
0.00 
0.00 
0.00 


-32.0416 

0.00 

-2.8771 

0.00 

0.0362 

0.00 

0.00 

0 

0 
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anc 


5  = 


-4.3850 

0 

0 

8.7745 

0.00 

5.6803 

0 

0 

-37.8929 

0 

0 

0 

0.00 

0.6216 

45.9851 

0 

-21.4821 

0 

0.00 

0 

0.00 

-7.1281 

-6.1465 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

Complete  model  linearized  with  0o  =  0  instead  of  0o  =  Qo  The  initial  conditions 


are  now: 


xO  = 


73.3000 
0 
0 
0 
0 
0 
0 
0 
0 


The  state  and  input  vectors 

obtained  from  trimming 

I  at  this  ! 

5tat 

e  are 

•  73.3000  ■ 

-0.00 

X  = 

1.6086 

-0.00 

-0.00 

0.00 

0.00 

ar 

id  u  = 

-0.0181  " 
-0.00 
0.00 
0.2336  . 

-0.00 

-0.00  . 

with  the  linearized  A 

and  B  matrices:  A  = 

z 

•  -0.0635 

0.00 

0.3277                0 

-1.4922 

0 

-0.00 

-32.1740    0 

-0.00 

-0.3911 

-0.00       1.6086 

0 

-72.6109 

32.1740 

-0.00    0 

-0.7572 

-0.00 

-4.7741                0 

67.9934 

0 

-0.0002 

-0.0002    0 

0.00 

-0.1471 

-0.00    -5.4414 

-0.00 

1.5183 

0.00 

-0.00    0 

0.0151 

-0.00 

-0.1933                0 

-3.1672 

0 

0.00 

-0.00    0 

0.00 

0.1440 

0.00    -1.0578 

0.00 

-0.8114 

0.00 

-0.00    0 

0 

0 

0           1.00 

0 

-0.00 

0 

-0.00    0 

0 

0 

0               0 

1.00 

-0.00 

-0.00 

0    0 

0 

0 

0 

0 

0.0 

0 

1.0 

0 

-0.00 

-0.00    0 

93 


B  = 


-4.5835 

0 

0 

8.7745 

0.00 

5.8282 

0 

0 

-38.8681 

0 

0 

0 

-0.00 

0.6252 

47.1717 

0 

-22.0417 

0 

0 

0 

-0.00 

-7.3151 

-6.4345 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

The  eigenvalues  for  the  case  where  ©o  =  0  are 

0 

-0.5285  -h  3.63462 

-0.5285  -  3.63462 

-5.6291 

eig{A)=      -3.9833  +  3.5521? 

-3.9833  -  3.5521? 

0.0420 

-0.0191  +  0.4963? 

-0.0191  -  0.4963? 


C.     CESSNA  172  RESULTS 

The  Cessnal  172  model  was  trimmed  and  linearized  using  the  state  vector 

219 
0 
0 
0 
2-0  =  0 

0 
0 
0 
0 

as  specified  in  [Ref.  Ros  79].  The  states  u,  u',  and  0  were  held  constant  making  the 
term  ?x  =  [1  3  8]'. 

Only  the  complete  model  was  linearized,  with  the  results  of  the  trim  routine 
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given  as: 

219.00 

0.00 

-0.0394 

-0.00 

X  =  —0.00       and  u 

-0.00 

-0.00 

-0.00 

0.00 

The  linearized  model  had  the  following  results  for  the  A  and  B  matrices  A  = 


0.0001 
0.00 
0.00 

1.0013 


and 


0.0442 

0.00 

0.0848 

0 

0.0382 

0 

0.00 

-32.1740 

0 

0.00 

-0.1620 

-0.00   - 

-0.0394 

0 

— 

217.2141 

32.1740 

-0.00 

0 

0.2916 

-0.00  - 

-2.1805 

0 

212.5399 

0 

-0.0002 

-0.0002 

0 

0.00 

-0.1313 

0.00  - 

12.4093 

0.00 

2.5342 

-0.00 

-0.00 

0 

0.0024 

-0.00  - 

-0.1085 

0 

— 

6.0778 

0 

0.00 

-0.00 

0 

-0.00 

0.0462 

0.00   - 

-0.3807 

0.00 

-1.2600 

0.00 

-0.00 

0 

0 

0 

0 

1.00 

0 

-0.00 

0 

-0.00 

0 

0 

0 

0 

0 

1.00 

-0.00 

0.00 

0 

0 

0 

0 

0 
■  -6.2509 

0 

0 

0.00 

0 

1.00 
3.2255  ■ 

-0.00 

-0.00 

0 

-0.00 

19.4571 

0 

0 

-44.3392 

0 

0 

0 

0.00 

4.7446 

57.4954 

0 

B  = 

-39.4919 

0 

0 

0 

5 

-0.00 

-10.2288 

-8.2562 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0  . 

and  the  following  eigenvalues 


eig{A)  = 


0 
-12.4309 
-0.6947  -f  3.30802 
-0.6947  -  3.3080? 
-4.1303 -h  4. 3895? 
-4.1303-4.3895? 
-0.0109 
-0.0209  +  0.1794? 
-0.0209-0.1794? 
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APPENDIX  C:  PROGRAM  LISTINGS 

A.     AROD  Matlab  ROUTINES 
1.     Main  Routine 

function  accel=main(vstate,m,rho,A,R) 

%  Function  will  compute  the  accelerations  due  to  the 

5i  gravitational  forces,  aerodynamic  forces  and  moments, 

%  and  control  forces  amd  moments. 

%  The  values  for  S,rho,m,b,8Lnd  c  are  used  as  input 

X  arguments  to  the  function  call,  and  are  loaded 

X  from  the  workspace.  There  should  be  a  file, 

%  loaddata.m  loaded  prior  to  running  the 

*/,  simulation. 

y,  define  states  in  terms  of  the  input  vector 

u=vstate(l) ; 

v=vstate(2) ; 

w=vstate(3) ; 

p=vstate(4) ; 

q=vstate(5) ; 

r=vstate(6) ; 

phi=vstate(7) ; 

theta=vstate(8) ; 

psi=vstate(9) ; 
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*/»     Define  the  control  inputs 
de=vstate(10) ; 
dr=vstate(ll) ; 
da=vstate(12) ; 
drpm=vstate(13) ; 


*/  this  subroutine  computes  linear  and  eingular  accelerations 

*/,  given  angular  and  linear  velocities; 

*/,  the  input  is  6x1  vector  =  [u  v  w  p  q  r]  ' 

y,  the  output  is  feedback  part  of  d/dt  [u  v  w  p  q  r]  ' 

5i  the  output  also  includes  the  Euler  angle  derivatives,  based  on 

y,  a  2-3-1  trernsf ormation,  for  Ldot,  used  m  the  function  grav.m. 

V  =  vstated  :3) ; 

omega  =  vstate(4:6); 

vdot  =  -crpr (omega, v) ; 

[Ibjir]  =  inertia; 

y,  compute  the  angular  momentum  due  to  the  body's  inertia,  lb 

Lb  =  lb  *omega; 

y,  compute  the  angular  momentum  due  the  spinning  prop's  inertia,  Ir 

0megaR=drpm*2*pi/60;  y,  angular  velocity  of  the  prop,  rad/sec 

Lr=Ir*[OmegaR;0;0]  ; 

temp=Lr+Lb; 

omdot  =  -Ib\crpr(omega,temp) ; 
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vstatedot=[vdot ;omdot] ; 

'/,  Use  the  Euler  angle  propagation  equations  for  a  2-3-1  rotation  sequence 
Ldot=[p/cos(psi)-cos(phi)*sin(psi)/cos(psi)*q+sin(phi)*sin(psi)/cos(psi)*r; 

cos(phi)/cos(psi)*q-sin(phi)/cos(psi)*r; 

sin(phi)*q+cos(phi)*r]  ; 


Given  the  vector  containing  the  state  derivatives, 

The  function  will  compute  the  forces  and  moments  due  to 

the  control  derivatives,  Cfd,  where  this 

is  dCf/dd. 

The  values  for  rho,A,R,m  are  used  as  input 

arguments  to  the  function  call,  and  are  loaded 

from  the  workspace.  There  should  be  a  file, 

arod.mat  loaded  prior  to  running  the 

simulation. 

hover  case  V=0,  dimensionalize  the  control  derivatives  based  on 
induced  velocity  through  the  rotor  disk,  Vi=sqrt (T/2/A/rho) 

Calculate  the  quantities  needed  for  the  force 

calculation. 

Call  a  function  to  return  the  stability  derivatives  wrt  to  moments 

Could  put  a  call  to  a  lookup  table  here 

Syntax— Z  =  TABLE2(TAB,X0.Y0)  or  Y  =  TABLE1(TAB,X0) 
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Cfd=getcfd; 

y.  Define  the  derivatives 

Clda=Cfd(4,3) 

Cmde=Cfd(5,l) 

Cndr=Cfd(6,2) 

y.  Calculate  the  Force  due  to  Cfd  derivatives 

y,  No  Aerodynamic  Forces  in  a  Hover 

D=0; 

Y=0 

L=0 

y,  calculate  the  force  due  to  thrust  in  body  coordinates 
y.  USING  THRUST  VS.  RPM  FROM  BOB  STONEY'S  TEST  RUNS 

T=0.0297*drpm-104.7; 

Vi=sqrt(T/2/A/rho) 

qbar= .5*rho*Vi"2;   y,  Vi  is  induced  velocity,  not  forward  speed 

Fout=[D;Y;L]; 

Fout = (Fout+ [T ; 0 ; 0] ) /m ; 

y.  Calculate  the  Moment  due  to  Cfd  derivatives 

y,  Itr  relates  the  duct  swirl  to  the  moment,  1,  produced  by  thrust 

ltr=-0.0542*T-0.9138; 

l=qbar*A*R* (Clda*da) +lt r ; 

m=qbar*A*R*(Cmde*de) ; 

n=qbar*A*R*(Cndr*dr) ; 
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Nout=Ib\[l;m;n] ; 
FNcfx=[Fout;Nout]  ; 


Given  the  vector  containing  the  state  derivatives, 

and  euler  angles,  the  function  will  compute 

the  forces  due  to  gravity  acting  on  the  aircraft. 


Calculate  gravitational  force,  based  on  a  2-3-1  Euler  angle 
rotation  for  position  of  the  aircraft.  Rotation  matrix  is  Ru2body 
Z  for  {U}  is  positive  down  and  X  for  {B}  is  straight  up. 


FNgrav=32.174*[-sin(theta)*cos(psi) ; 

sin(theta)*sin(psi)*cos(phi)+sin(phi)*cos(theta) ; 
cos(phi)*cos(theta)-sin(theta)*sin(psi)*sin(phi) ; 
0;0;0]; 

y,       Sura  up  the  normalized  forces  eind  moments  to  feed  back  into  the  integrator 

vstatedot= [vstatedot+FNcf x+FNgrav ; Ldot] ; 
accel=vstatedot ; 

2.      Supporting  Subroutines 

function  y  =  crpr(omega,x) 
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%   this  subroutine  computes  the  crossproduct  of  omega  and  x: 
•/,  y  =  omega  X  x 
p  =  omega(l) ; 
q  =  omega(2) ; 
r  =  omega(3) ; 
t  =  [0  -r  q; 
r  0  -p; 
-q  p  0]  ; 
y  =  t  *  x; 


function  Fgrav=grav(x) 

•/,  GRAY  will  compute  the  gravitational  force 

y,  acting  on  the  body,  in  body  coordinates 

g=x(l); 

phi=x(2); 

theta=x(3) ; 

psi=x(4) ; 

Fgrav=g*[-sin(theta)*cos(psi) ; 

sin(theta)*sin(psi)*cos(phi)+sin(phi)*cos(theta) ; 

cos(phi)*cos(theta)-sin(theta)*sin(psi)*sin(phi)] ; 


function  [ib,ir]  =  inertia 

*/,  this  subroutine  creates  inertia  matrices  called  ib 
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*/,  and  ir  for  the  body  and  rotor  inertia,  respectively, 

*/,  Arod  hover  case 

ixx  =  1.2312; 

iyy  =  3.9584; 

izz  =  3.9825; 

ixz  =  0; 

irx=. 00898; 

iry=.0045; 

irz=.0045; 

ib  =  [ixx  0  -ixz;0  iyy  0;-ixz  0  izz]; 

ir  =  [irx  0  0;0  iry  0;  0  0  irz]  ; 

3.  Data  and  Initialization  Subroutines 

'/,  load  bluebird  data 

A=3.14;  y.  Area  of  rotor  disk,  ft"2 

Vt=712.09;  '/,  Rotor  tip  speed,  rad/sec 

m=2.6419;  •/.  Mass,  slugs 

R=l;  y,  Radius  of  Rotor  Blade,  ft 

rho=. 002377;  */,   Air  density,  slug/ft "3 

Uo=0;  y.  Nominal  Velocity,  ft/sec 

y.  Initial  Euler  Angle  Orientation,  radians 

Phi=0; 

Theta=1.57;  %  Same  as  Steady  State  alpha 

Psi=0; 

Lo=[Phi;Theta;Psi] ; 
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*h     Initial  Conditions  to  determine  the  Aircraft 

y.  Linear  and  Rotational  Velocity  States 

y,  u,v,w  are  computed  from  Uo,  alpha,  and  beta 

y,  p,q,r  are  assumed  as  zero. 

alpha=Theta; 

beta=0; 

Xo= [Uo ; alpha ; beta] ; 

y.  Returns  Initial  Conditions  for  the  main 
y,  integrator  in  the  rigid  body  EOM  block 
iO=init_var(Lo,Xo) ; 


function  Cfd=getcfd 

Cfd=getcfd  will  return  values  for 
The  stability  derivatives  for  D,Y,L,l,m,and  n 
due  to  the  control  inputs, 
format  for  data  is; 
[CDde  CDdr  CDda 
CYde  . . . 
CLde  . . . 
Clde  Cldr  Clda 
Cmde  . . . 
Cnde  . . .   Cnda 
C  Data  is  non-dimensionalized  by  using  induced  velocity 
C  Vi=sqrt(T/2/A/rho) 
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*/,  for  arod  hover  case 
Cfd=[    0    0    0; 
0    0    0; 
0    0    0 
0    0   1.438; 
-1.233  0     0; 
0  -1.233  0] ; 


function  iO=init_var(Lo,Xo) 

y,  INIT_VAR(X)  will  initialize  the  integrators 

*/,  initial  states,  iO,  given  the  initial 

%  conditions  desired. 

y.  Required  initial  conditions  are  the  Euler 

y,  ajigle  orientation,  total  velocity,  Uo,  initial 

y,  ADA,  cind  sideslip  angle,  beta. 

y,  Lo=[phi  ;  theta;psi]  ' 

y.  Xo=[Uo;  alpha;  beta]  ' 

y.  All  body  rotation  rates  are  assumed  to  be  zero 


y.  Initialize  the  states  u,v,w,p,q,r 

Uo=Xo(l); 

alpha=Xo(2) ; 

beta=Xo(3); 
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Ca=cos (alpha) ; 
Sa=siii (alpha)  ; 
CB=cos(beta) ; 
SB=sin(beta) ; 

Rwb=[Ca*CB  -Ca*SB  -Sa;SB  CB  0;Sa*CB  -Sa*SB  Ca] ; 

vO=Rwb*[Uo;0;0]; 

wO=[0;0;0]; 

iO=[vO;wO;Lo] ; 


function  Qo=initq(lainbda) 

y.  Function  initQ  will  return  values  for 

*l,     the  quaternion  DCM  based  on  a  given 

*/  set  of  Euler  atngles . 

y.  Set  for  a  Euler  3-2-1  rotation 

y,  C)(1)=B0 

y.  Q(2)=B1 

y.  Q(3)=B2 

y.  Q(4)=B3 

phi=lainbda(l) ; 

theta=leLmbda(2) ; 

psi=lambda(3) ; 

Qo(l)=.5*sqrt(l+cos(psi)*cos(theta)+sin(psi)*sin(theta)*sin(phi) 
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+cos(psi)*cos(phi)+cos(theta)*cos(phi)) ; 
Qo(2)=l/4/Qo(l)*(cos(theta)*sin(phi)-sin(psi)*sin(theta)*cos(phi) 

+cos(psi)*sin(phi)) ; 
Qo(3)=l/4/Qo(l)*(cos(psi)*sin(theta)*cos(phi)+sin(psi)*sin(phi) . . 

+sin(theta) ) ; 
Qo(4)=l/4/Qo(l)*(sin(psi)*cos(theta)-cos(psi)*sin(theta)*sin(phi) 

+sin(psi)*cos(phi)) ; 
Qo=Qo'  ; 


B.     BLUEBIRD  Matlab  ROUTINES 
1.     Main  Routine 

function  accel=iiiain(vstate,rho,b,c,S,m,Xo) 

Function  will  compute  the  accelerations  due  to  the 
gravitational  forces,  aerodynconic  forces  and  moments, 
and  control  forces  and  moments. 

The  values  for  S,rho,m,b,and  c  are  used  as  input 
arguments  to  the  function  call,  and  are  loaded 
from  the  workspace.  There  should  be  a  file, 
loaddata.m  loaded  prior  to  running  the 
simulation. 


y,  define  states  in  terms  of  the  input  vector 
u=vstate(l) ; 
v=vstate(2) ; 
w=vstate(3) ; 
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p=vstate(4) ; 
q=vstate(5)  ; 
r=vstate(6) ; 
phi=vstate(7)  ; 
theta=vstate(8) ; 
psi=vstate(9) ; 


X  Define  the  control  inputs 

de=vstate(10) 

dr=vstate(ll) 

da=vstate(12) 

dt=vstate(13) 


y,  calculate  the  aerodynamic  terms 
Vt=sqrt(u'2+v"2+w'2) ; 
qbar= .  5*rho*  (Vt)  "^2 ; 

Ib=inertia; 

*/,  wind  to  body  transformation 
Rwb=rw2b(vstate,Vt) ; 


y. 
y. 
y.- 
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*/,  CHI  will  compute  the  left  hajid  side  of  the  state 

y,  equation.  This  is  the  term  dependant  on  dCf/dxdot. 

y.  Now  calculate  the  S  matrix  that  non-dimensionalizes 

y,  the  moments.  Also  includes  the  correction  for  Lift  and  Drag 

y,  acting  in  the  direction  opposite  to  the  positive  coordinate 

y,  direction. 

*/t   Get  the  stability  derivatives  for  forces  and  moments 

Cf xdot=getcf xd ; 

CLad=Cfxdot(3) ; 

Cmad=Cfxdot(5) ; 

y,  Twb  is  a  intermediate  step 

Twb=[eye(3)/m  zeros(3) ;zeros(3)  inv(Ib)]*[Rwb  zeros (3) ; zeros (3)  Rwb] ; 

y,  calculate  the  M2  matrix  to  allow  use  of  the  states,  rather  than 
y,  the  normalized  states.  Only  accounts  for  the  alpha  dot  variables 
y,  since  the  beta  dot  terms  are  not  ordinarily  computed. 

y.M2=[0  0  10  0  0]*c/2/Vt^2; 

y,Sprime=diag([-S  S  -S  S*b  S*c  S*b]); 

y.  To  save  some  math  here,  the  product  of  Sprime,  Cfxdot,  and  M2  is: 

PROD=[0  00000;000000;00  -S*CLad*c/2/Vt-2  0  0  0; 

0  0  0  0  0  0;0  0  S*c*Cmad*c/2/Vt^2  0  0  0;0  0  0  0  0  0]; 
y.  Calculate  chi 
chil=(eye(6)-Twb*qbar*PR0D) ; 

y. 

108 


Given  the  vector  containing  the  state  derivatives, 

sind  euler  angles,  the  function  will  compute 

the  forces  due  to  gravity  acting  on  the  aircraft. 

Calculate  gravitational  force,  based  on  a  3-2-1  Euler  angle 
rotation  for  position  of  the  aircraft.  Rotation  matrix  is  Ru2body 


Fgrav=32 . 174* [-sin(theta) ; 

cos(theta)*sin(phi) ; 
cos(theta)*cos(phi)] ; 

y,  Premultiply  by  the  Chi~-1  term  from  the  left  hand  side 
FNgrav=chi 1\ [Fgrav ; 0 ; 0 ; 0] ; 

y. 
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y,  Cfx(u)  Given  the  vector  containing  the  state  derivatives, 

y,  The  function  will  compute  the  forces  and  moments  due  to 

y,  the  stability  derivatives,  Cfx',  where  this 

y,  is  dCf/dx. 

y. 
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*/,  Call  a  function  to  return  the  stability  derivatives  wrt  to  moments 

*/,  Could  put  a  call  to  a  lookup  table  here 

•/.  Syntax— Z  =  TABLE2(TAB,X0,Y0)  or  Y  =  TABLE1(TAB,X0) 

Cfx=getcfx; 

CDu=Cfx(l,l);  CDa=Cfx(l,3); 

CYb=Cfx(2,2);  CYp=Cfx(2,4) ;  CYr=Cfx(2,6) ; 

CLu=Cfx(3,l);  CLa=Cfx(3,3);  CLq=Cfx(3,5) ; 

Clb=Cfx(4,2);  Clp=Cfx(4,4) ;  Clr=Cfx(4,6) ; 

Cmu=Cfx(5,l) ;  Cma=Cfx(5,3) ;  Cmq=Cfx(5,5) ; 

Cnb=Cfx(6,2);  Cnp=Cfx(6,4) ;  Cnr=Cfx(6,6) ; 


ss=getcf 0 
CDO=ss(l) 
CL0=ss(3) 
Cm0=ss(5) 


Cfd=getcfd; 

*/,  Define  the  derivatives 

CDde=Cfd(l,l); 

CYdr=Cfd(2,2);  CYda=Cf d(2,3) ; 

CLde=Cfd(3,l); 

Cldr=Cfd(4,2);  Clda=Cf d(4,3) ; 

Cmde=Cfd(5,l); 

Cndr=Cfd(6,2);  Cnda=Cf d(6,3) ; 
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y.  Calculate  the  Force  due  to  Cfx'  derivatives 
*/,     And  the  control  derivatives 
*/,  in  the  wind  coordinates 

D=-S*qbar/m* (CDO+CDa*w/Vt+CDde*de) ; 

Y=S*qbar/m* (CYb*v/Vt+CYr*r*b/2/Vt+CYdr*dr+CYda*da) ; 

L=-S*qbar/m* (CL0+CLa*w/Vt+CLq*q*c/2/Vt+CLde*de) ; 

5i  calculate  the  force  due  to  thrust  in  body  coordinates 
y.  THRUST  IS  ESTIMATED,  BASED  ON  4.0  HP,  PROP  EFFICENCY=.65 

T=15; 
Xt=T/m*dt ; 

Fout=Rwb*([D;Y;L]); 
Fout=Fout+[Xt;0;0] ; 

y.  Calculate  the  Moment  due  to  Cfx'  derivatives 
y.  And  the  control  derivatives 

l=qbar*S*b*(Clb/Vt*v+Clp*b/2/Vt*p+Clr*b/2/Vt*r+Cldr*dr+Clda*da) ; 
m=qbar*S*c* (Cm0+Cma/Vt*w+Cmq/Vt*c/2*q+Cmde*de) ; 
n=qbar*S*b*(Cnb/Vt*v+Cnp*b/2/Vt*p+Cnr*b/2/Vt*r+Cndr*dr+Cnda*da) ; 

Nout=Ib\ (Rwb* [1 ; m ; n] ) ; 
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%  Premultiply  by  the  Chi"-1  term  from  the  left  hand  side 

FNcfx=chil\[Fout;No\it]  ; 

% 

% 


%   this  subroutine  computes  linear  ajid  angular  accelerations 

%   given  emgular  cind  linear  velocities; 

*/,  the  input  is  6x1  vector  =  [u  v  w  p  q  r]  ' 

y,  the  output  is  feedback  part  of  d/dt  [u  v  w  p  q  r]  ' 

%  the  output  also  includes  the  Euler  angle  derivatives 

*/,  Ldot,  used  in  the  function  grav.m. 

v  =  vstated  :3) ; 

omega  =  vstate(4:6); 

vdot  =  -crpr (omega, v) ; 

temp  =  Ib*omega; 

omdot  =  -Ib\crpr(omega,temp) ; 

vstatedot=chil\[vdot ;omdot]  ; 

•/,  Use  the  Euler  angle  propagation  equations 
Ldot=[p+sin(phi)*tan(theta)*q+cos(phi)*tcLn(theta)*r; 

cos(phi)*q-sin(phi)*r; 

sin(phi)/cos(theta)*q+cos(phi)/cos(theta)*r] ; 
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vstatedot= [vstatedot+FNcf x+FNgrav ; Ldot] ; 
accel=vstatedot ; 


2.      Supporting  Subroutines 

function  y  =  crpr(omega,x) 

%   this  subroutine  computes  the  crossproduct  of  omega  and  x: 
y,  y  =  omega  X  x 
p  =  omega(l) ; 
q  =  omega(2) ; 
r  =  omega(3) ; 
t  =  [0  -r  q; 
r  0  -p; 
-q  p  0]  ; 
y  =  t  *  x; 


function  Fgrav=grav(x) 

y,  GRAY  will  compute  the  gravitational  force 

y.  acting  on  the  body,  in  body  coordinates 

g=x(l); 
phi=x(2); 
theta=x(3) ; 
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Fgrav=g* [-sin(theta) ; 

cos(theta)*sin(phi) ; 
cos(theta)*cos(phi)]  ; 


function  lb  =  inertia 

y,  this  subroutine  creates  inertia  matrix  called  lb 

*/,  for  the  Bluebird  test  aircraft. 

'/,  All  units  are  slug-ft"2 

Ix=10.0; 

Iy=16.12; 

Iz=7.97; 

Ib=[Ix  0  0;0  ly  0;0  0  Iz]  ; 


function  [Rwb]=Rw2b(x,Vt) 

*/,  RWIND2BQDY  Rotation  matrix  for  wind  to  body  coordinate  trainsf ormations 

Ca=x(l)/Vt; 
Sa=x(3)/Vt; 
CB=x(l)/Vt; 
SB=x(2)/Vt; 

Rwb=[Ca*CB  -Ca*SB  -Sa;SB  CB  0;Sa*CB  -Sa*SB  Ca] ; 
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3.     Data  and  Initialization  Subroutines 

y,  load  bluebird  data 


S=22.38;  " 

Uo=73.3; 

m=1.7095; 

b=12.42; 

c=1.802; 

rho=. 002377; 


Planform  Area,  ft"2 

Nominal  Velocity,  ft/sec 

Mass,  slugs 

Span,  ft 

Average  Chord,  ft 

Air  density,  slug/ft'3 


y,  Initial  Euler  Angle  Orientation,  radians 

Phi=0; 

Theta=.0912;         y.  Same  as  Steady  State  alpha 

Psi=0; 

Lo=[Phi;Theta;Psi]  ; 

y.  Initial  Conditions  to  determine  the  Aircraft 

y,  Linear  and  Rotational  Velocity  States 

%     u,v,w  are  computed  from  Uo,  alpha,  emd  beta 

y,  p,q,r  are  assumed  as  zero. 

alpha=Theta; 

beta=0; 

Xo=[Uo;alpha;beta]  ; 

y.  Returns  Initial  Conditions  for  the  main 
y,  integrator  in  the  rigid  body  EOM  block 
iO=init_var(Lo,Xo) 
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function  iO=init_var(Lo,Xo) 

y,  INIT_VAR(X)  will  initialize  the  integrators 

*/,  initial  states,  iO,  given  the  initial 

7,  conditions  desired. 

X  Required  initial  conditions  are  the  Euler 

*/,  angle  orientation,  total  velocity,  Uo,  initial 

*/,  ADA,  and  sideslip  angle,  beta. 

*/,  Lo=[phi  ;theta;psi]  ' 

*/.  Xo=[Uo;alpha;beta]  ' 

y.  All  body  rotation  rates  are  assumed  to  be  zero 

y.  Initialize  the  Euler  eingle  DCM 

y.  Initialize  the  states  u,v,w,p,q,r 
Uo=Xo(l); 
alpha=Xo(2) ; 
beta=Xo(3) ; 

Ca=cos (alpha) ;  , 

Sa=sin(alpha) ; 
CB=cos(beta)  ; 
SB=sin(beta)  ; 

Rwb=[Ca*CB  -Ca*SB  -SaiSB  CB  0;Sa*CB  -Sa*SB  Ca]  ; 
vO=Rwb*[Uo;0;0]  ; 
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wO=[0;0;0]; 
iO=[vO;wO;Lo] ; 


function  CfO=getcfO 

*h   CfO=getcfO  will  return  values  for 

y,  the  nominal  values  for  coefficients 

*/,  format  of  input  is  [CDO  CYO  CLO  CIO  CmO  CnO]  '  ; 

Cf0=[0.03  0  0.3  0  0  0]  '; 


function  Cfd=getcfd 

Cfd=getCfd_F(n)  will  return  values  for 
The  stability  derivatives  for  D,Y,cmd  L 
due  to  the  control  inputs, 
format  for  data  is; 
[CDde  CDdr  CDda 

CYde  . . . 

CLde  . . . 

Clde  Cldr  Clda 

Cmde  . . . 

Cnde  .  .  .  Cnda] 

For  the  test  aircraft  Bluebird 
C  Derivatives  from  DATCOM 
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Cfd=[   .065   0     0; 
0    .0697   0; 
.472   0     0 

0   .0028  .265; 

-1.41  0  0; 

0  -.0329  -.0347]; 


function  Cfx=getcfx 

*/,  Cfx=getcfx(n)  will  return  values  for 

*/,  The  stability  derivatives  for  D,Y,and  L 

y,  due  to  the  state  vector. 

y,  format  of  data: 

y.[CDu  CDb  CDa  CDp  CDq  CDr; 

y.  CYu  ... 

y.  cLu  ... 
y.  ciu  ... 

y,  Cmu  .  .  . 

y.  Cnu  .  .  .  ] 

y.  For  the  test  aircraft  Bluebird 
y.  Derivatives  from  DATCOM 

Cfx=[0    0   .188     0    0     0; 
0  -0.31   0      0    0   .0973; 
0    0    4.22    0   3.94    0; 
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0  -.0597  0  -.363  0  .1; 
0  0  -1.163  0  -11.77  0; 
0   .0487   0  -.0481  0    -.0452]; 


function  Cfxdot=getcfxd 

Cfxdot=getcfxd(n)  will  return  values  for 
The  stability  derivatives  for  D,Y,aiid  L 
due  to  the  state  vector.  Beta  dot  terms  are  ignored 
since  they  are  not  normally  determined, 
format  is: 
[  CDadot 

CYad 

CLad 

Clad 

Cmad 

Cnad  ] 


Cf xdot=  [0 ; 0 ; 1 . 32 ; 0 ; -4 . 7 ; 0]  ; 


119 


REFERENCES 

[Siu  91]  Siuru,  W.D.,  Planes  Without  Pilots:  Advances  in  Unmanned  Flight,  TAB 
Books,  Blue  Ridge  Summit,  PA,  1991. 

[Wh  87]  White,  J.E.,  and  Phelan,  J.R.,  "Stability  Augmentation  for  a  Free  Flying 
Ducted  Fan,"  Proceedings  of  the  AIAA  Guidance,  Navigation,  and  Control  Con- 
ference, Monterey,  CA,  Aug.  1987,  pp  896-904. 

[MCG  87]  Jennings  Jr,  J.F.,  "Why  the  Corps  Needs  Robots,"  Marine  Corps  Gazette, 
Vol.  5,  No.  1987,  pp  36-39. 

[Sa  89]  Not  Attributed,  Sandia  Science  News  Vol.  1,  1989 

[Kre  92]  Kress,  G.A.,  "Preliminary  Development  of  a  VTOL  Unmanned  Air  Vehicle 
for  the  Close-Range  Mission."  Master's  Thesis,  Department  of  Aeronautics,  Naval 
Postgraduate  School,  Monterey,  CA,  1992. 

[Sto  93]  Stoney,  R.B.,  "Design,  Fabrication,  and  Test  of  a  Vertical  Attitude  Take-Off 
and  Landing  Unmanned  Air  Vehicle,"  Engineer's  Thesis,  Department  of  Aero- 
nautics, Naval  Postgraduate  School,  Monterey,  CA,  June  1993. 

[We  88]  Weir,  R.J.,  "Aerodynamic  Design  Considerations  for  Free-Flying  Ducted 
Propellor,"  Proceedings  of  the  1988  Atmospheric  Flight  Mechanics  Conference, 
AIAA,  Washington,  D.C.,  Aug.  1988,  pp.  720-731. 

[DOD  92]  "DoD  Unmanned  Aerial  Vehicle  Master  Plan,"  Department  of  Defense, 
Washington,  D.C.,  1992. 

[Wh  91]  White,  J.E.,  and  Phelan,  J.R.,  "Stability  Augmentation  and  Control  Decou- 
pling for  the  Airborne  Remotely  Operated  Device,"  Journal  of  Guidance,  Control, 
and  Dynamics,  Vol.  14,  No.l,  1991,  pp  176-183. 

[Sil  91]  Silvestre,  C.J.,  "Modeling  and  Control  of  Underwater  Vehicles,"  Master's 
Thesis,  Department  of  Electrical  Engineering,  Institute  Superior  Tecnico,  Lisbon, 
Portugal,  1991. 

[Cra  86]  Craig,  J.J.,  Introduction  to  Robotics  Mechanics  and  Control  Addison- 
Wesley,  New  York,  1986. 

[Ju  92]  Junkins,  J.L.,  An  Introduction  to  Dynamics  and  Control  of  Fleiihle  Struc- 
tures, AIAA,  Washington  D.C.,  1992. 

[Sch  92]  Schmidt,  L.V.,  Class  Notes  for  AE3340,  U.S.  Naval  Postgraduate  School, 
Monterey,  CA.  1992. 

[Ka  83]  Kane,  T.R.,  Likins,  P.W.,  Levinson,  D.A.,  Spacecraft  Dynamics,  McGraw- 
Hill,  New  York,  1983. 


120 


[Mo  84]  Morton,  H.S.,  "A  Formulation  of  Rigid-Body  Rotational  Dynamics  Based  on 
Generalized  Angular  Momentum  Variables  Corresponding  to  the  Euler  Parame- 
ters," Proceedings  of  the  AIAA/AAS,  Seattle,  WA,  August  1984. 

[Ro  58]  Robinson,  A.C.,  "On  the  Use  of  Quaternions  in  Simulation  of  Rigid  Body 
Motion,"  WADC  Technical  Report  58-17,  Wright  Air  Development  Center,  De- 
cember, 1958. 

[Whi  59]  Whittaker,  E.T.,  A  Treatise  on  the  Analytical  Dynamics  of  Particles  and 
Rigid  Bodies,  Cambridge  Univ.  Press,  4th  Edition,  1959. 

[Gre  88]  Greenwood,  D.T.,  Principles  of  Dynamics,  2nd  Ed.,  Prentice-Hall,  Engle- 
wood  Cliffs,  N.J.,  1988. 

[Ros  79]  Roskam,  J.,  Airplane  Flight  Dynamics  and  Automatic  Flight  Controls, 
Roskam  Aviation  and  Engineering  corp,  Ottawa,  KS,  1979 

[Th  89]  Thompson,  CM.,  "Aircraft  Equations  of  Motion  and  Forming  Linear  Mod- 
els," Boeing  Document  D6-54972,  Boeing  Commercial  Airplane  Company,  Seat- 
tle, Washington,  1989. 

[USAF  60]  USAF  STABILITY  AND  CONTROL  HANDBOOK,  Wright  Air  Devel- 
opment Division,  United  States  Air  Force,  Wright  Patterson  AFB  McGregor  and 
Werner,  Inc.,  Dayton,  OH,  1960. 

[Pro  90]  Prouty,  R.W.,  Helocopter  Performance,  Stability,  and  Control,  Robert  E. 
Krieger,  Malabar,  Florida,  1990. 

[WAT  93]  Watson  Industries,  Techical  Specifications  for  IMU-600D  Watson  Indus- 
tries, Eau-Claire,  WI,  1993. 

[Bro  64]  Broxmeyer,  C,  Inertial  Navigation  Systems,  McGraw-Hill,  New  York,  1964. 

[St  88]  Strum,  R.D.,  and  Kirk,  D.E.,  First  Principles  of  Discrete  Systems  and  Digigal 
Signal  Processing,  Addison-Wesley,  New  York,  1988. 

[Sp  89]   Spong,  Vydia,  Sagas.,  Robot  Dynamics  and  Control,  Wiley,  New  York,  1989. 


121 


INITIAL  DISTRIBUTION  LIST 


No.  of  Copies 

1.  Defense  Technical  Information  Center  2 
Cameron  Station 

Alexandria,  Virginia  22304-6145 

2.  Commandant  of  the  Marine  Corps  1 
Code  TE06 

Headquarters,  U.S.  Marine  Corps 
Washington,  D.C.  20380-0001 

3.  Library,  Code  52  2 
Naval  Postgraduate  School 

Monterey,  California  93943-5002 

4.  Dr.  Isaac  I.  Kaminer  5 
Department  of  Aeronautics  and 

Astronautics,  Code  AA/Ka 
Naval  Postgraduate  School 
Monterey,  California  93943-5000 

5.  Dr.  Richard  W.  Howard  2 
Department  of  Aeronautics  and 

Astronautics,  Code  AA/Ho 
Naval  Postgraduate  School 
Monterey,  California  93943-5000 

6.  Chairman  2 
Department  of  Aeronautics  and 

Astronautics 

Naval  Postgraduate  School 

Monterey,  California  93943-5000 

7.  Capt.  David  R.  Kuechenmeister  2 
1995  Skidmore  Circle 

Lawrenceville,  Georgia  30244 


tf^P^  '/^ 


122 


DUDLEY  KNOX  LIBRARY 
NAVAL  POSTGRADUATE  SCHOOI 
MONTEREY  CA  93943-5101 


RES 

STSC?  suide lines  for 
\STSC, 

route  to:     RESERlJE 

Hahr?  Randoiph  L. 

10:32763000191530 


GAYLORD  S 


W: 


